{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 作业内容：\n",
    "\n",
    "> 编程完成3阶辛普森积分与复化辛普森积分，以如下函数为例，与sympy的结果对比\n",
    ">> $\\int_{-1}^1 \\sqrt{x - 1} dx$\n",
    "\n",
    "> 计算理想强子共振气体的某些物理量，并与论文数据对比"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 参考答案"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 1， 辛普森积分"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "三阶辛普森积分公式如下：\n",
    "\n",
    "$$F = (b - a)\\left[ {1\\over 8} f(a) + {3 \\over 8} f({2a + b \\over 3}) + {3 \\over 8} f({a + 2 b \\over 3}) + {1 \\over 8} f(b) \\right]$$\n",
    "\n",
    "其中，b与a分别为积分上下限。\n",
    "\n",
    "那么显然可得："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simpson result: 2.000, analytical solutions: 2.000\n"
     ]
    }
   ],
   "source": [
    "def MySimpson(f, a, b):\n",
    "    return (b - a) * (f(a) + 3*f(2*a + b) + 3*f(2*b + a) + f(b)) / 8\n",
    "\n",
    "#test using a simple function\n",
    "def demo_func(x):\n",
    "    return x**2\n",
    "print('Simpson result: {:.3f}, analytical solutions: {:.3f}'.format(MySimpson(demo_func, 0, 1), 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sympy as sym\n",
    "import cmath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x):\n",
    "    return (x - 1)**0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sympy result\n",
    "x, a, b = sym.symbols('x a b')\n",
    "sym_res = complex(sym.integrate(func(x), (x, -1, 1)))\n",
    "#My Simpson result\n",
    "a = -1 + 0j\n",
    "b = 1 + 0j\n",
    "my_res = MySimpson(func, a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Translate into polar cordinates\n",
    "sym_polar_res = cmath.polar(sym_res)\n",
    "my_polar_res = cmath.polar(my_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "复化辛普森积分就是把积分区间划分成n个子区间然后分别做辛普森积分，那么显然："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MyMultiSimpson(f, a, b, n, SimpsonMethod):\n",
    "    nodes = np.linspace(a, b, n)\n",
    "    nodes_left = nodes[:-1]\n",
    "    nodes_right = nodes[1:]\n",
    "    res = [SimpsonMethod(f, item_left, item_right) for item_left, item_right in zip(nodes_left, nodes_right)]\n",
    "    return np.sum(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_multi_res = MyMultiSimpson(func, a, b, 1000, MySimpson)\n",
    "my_multi_polar_res = cmath.polar(my_multi_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian\n",
      "Sympy result: 1.8856180831641267j\n",
      "My result: (8.659560562354934e-17+1.4142135623730951j)\n",
      "My Multi-Simpson result: (0.4714095442741609+1.8047382776107321j)\n",
      "Polar\n",
      "Sympy result: (1.8856180831641267, 1.5707963267948966)\n",
      "My result: (1.4142135623730951, 1.5707963267948966)\n",
      "My Multi-Simpson result: (1.8652901139249423, 1.3152984048101044)\n"
     ]
    }
   ],
   "source": [
    "#Comparison\n",
    "print('Cartesian\\nSympy result: {}\\nMy result: {}\\nMy Multi-Simpson result: {}'.format(sym_res, my_res, my_multi_res))\n",
    "print('Polar\\nSympy result: {}\\nMy result: {}\\nMy Multi-Simpson result: {}'.format(sym_polar_res, my_polar_res, my_multi_polar_res))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 2，状态方程"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "前面大部分内容是不用改动的，我这里直接把有用的复制过来了"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Unchanged part\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "from scipy.integrate import quad\n",
    "\n",
    "def dist_f(k, T, hadron):\n",
    "    '''the distribution function of hadron \n",
    "    :k: float or np.array(), momentum in units [GeV]\n",
    "    :T: float, temperature, [GeV]\n",
    "    :hadron: dict with {'mass', '2*spin+1', 'eta',\n",
    "                        'baryon_num', 'strange_num', 'charge_num'}\n",
    "             where eta: +1 for fermion, -1 for boson\n",
    "             \n",
    "    :return: (2s + 1) / (np.exp(np.sqrt(k**2 + m**2)/T) + eta)\n",
    "    '''\n",
    "    spin_dof = hadron['2*spin+1']\n",
    "    m = hadron['mass']\n",
    "    eta = hadron['eta']\n",
    "    return  spin_dof / (np.exp(np.sqrt(k**2 + m**2)/T) + eta)\n",
    "\n",
    "def density_n(T, hadron):\n",
    "    '''calc density for given temperatures\n",
    "    :T: float, temperature, [GeV]\n",
    "    :hadron: dict with {'mass', '2*spin+1', 'eta',\n",
    "                        'baryon_num', 'strange_num', 'charge_num'}\n",
    "             where eta: +1 for fermion, -1 for boson\n",
    "    :return: hadron density at T '''\n",
    "    # 此处 intg_n 是 inline function\n",
    "    intg_n = lambda k: k**2 * dist_f(k, T, hadron) / (2 * np.pi**2)\n",
    "    \n",
    "    kmax = 50 * T\n",
    "    # 不要直接使用 quad(f, 0, 50*T), 给 50*T 命名为 kmax，程序可读性更好\n",
    "    ndensity = quad(intg_n, 0, kmax)[0]\n",
    "    \n",
    "    return ndensity\n",
    "\n",
    "def pressure_p(T, hadron):\n",
    "    '''HRG presure for given temperatures\n",
    "    :T: float or np.array(), temperature, [GeV]\n",
    "    :hadron: dict with {'mass', '2*spin+1', 'eta',\n",
    "                        'baryon_num', 'strange_num', 'charge_num'}\n",
    "             where eta: +1 for fermion, -1 for boson\n",
    "    :return: hadron presure at T  in unit [GeV]^4 '''\n",
    "    spin_dof = hadron['2*spin+1']\n",
    "    m = hadron['mass']\n",
    "    eta = hadron['eta']\n",
    "    intg_p = lambda k:  k**2 * np.log(1 + eta * np.exp(-(np.sqrt(k**2 + m**2))/T))    \n",
    "    coef_p = eta * T * spin_dof / (2 * np.pi**2)\n",
    "\n",
    "    kmax = 50 * T\n",
    "    pressure = coef_p * quad(intg_p, 0, kmax)[0] \n",
    "    return pressure\n",
    "\n",
    "def energy_density_e(T, hadron):\n",
    "    '''HRG presure for given temperatures\n",
    "    :T: float or np.array(), temperature, [GeV]\n",
    "    :hadron: dict with {'mass', '2*spin+1', 'eta',\n",
    "                        'baryon_num', 'strange_num', 'charge_num'}\n",
    "             where eta: +1 for fermion, -1 for boson \n",
    "    :return: hadron energy density at T in unit [GeV]^4'''\n",
    "    m = hadron['mass']\n",
    "\n",
    "    intg_e = lambda k: k**2 * np.sqrt(k**2 + m**2) * dist_f(k, T, hadron) / (2 * np.pi**2)\n",
    "    kmax = 50 * T\n",
    "    edensity = quad(intg_e, 0, kmax)[0] \n",
    "    return edensity\n",
    "\n",
    "def eos(T, mu, hadron):\n",
    "    '''calc the pressure vs energy density\n",
    "    :T: float, temperature, [GeV]\n",
    "    :mu: (mu_B, mu_S, mu_Q)\n",
    "    :hadron: dict with {'mass', '2*spin+1', 'eta',\n",
    "                        'baryon_num', 'strange_num', 'charge_num'}\n",
    "             where eta: +1 for fermion, -1 for boson\n",
    "    return: pressure, energy density in unit [GeV]^4，particle density'''\n",
    "    pressure = pressure_p(T, hadron)\n",
    "    edensity = energy_density_e(T, hadron)\n",
    "    ndensity = density_n(T, hadron)\n",
    "    \n",
    "    return pressure, edensity, ndensity\n",
    "\n",
    "pdg_ = []\n",
    "name_ = []\n",
    "mass_ = []\n",
    "decay_width_ = []       # decay width\n",
    "spin_degeneracy_ = []    # 2*s + 1\n",
    "isbaryon_ = []\n",
    "strange_num_ = []\n",
    "charge_num_ = []\n",
    "with open(\"data/pdg05.dat\", \"r\") as pdg:\n",
    "    lines = pdg.readlines()\n",
    "    row = 0\n",
    "    while row < len(lines):\n",
    "        # 按空格将不同的数据分开\n",
    "        particle = lines[row].split()\n",
    "        # 最后一列存储了衰变通道个数\n",
    "        decay_channels = int(particle[-1])\n",
    "        row += decay_channels\n",
    "        \n",
    "        #print(particle[1], particle[10])\n",
    "        pdg_.append(int(particle[0]))\n",
    "        name_.append(particle[1])\n",
    "        mass_.append(float(particle[2]))\n",
    "        decay_width_.append(float(particle[3]))\n",
    "        spin_degeneracy_.append(int(particle[4]))\n",
    "        isbaryon_.append(int(particle[5]))\n",
    "        strange_num_.append(int(particle[6]))\n",
    "        charge_num_.append(int(particle[10]))\n",
    "        \n",
    "        row += 1\n",
    "        \n",
    "df = pd.DataFrame({\"pdg\":pdg_, \n",
    "                   \"name\":name_, \n",
    "                   \"mass\":mass_,\n",
    "                   \"decay_width\":decay_width_,\n",
    "                   \"2*spin+1\":spin_degeneracy_,\n",
    "                   \"isbaryon\":isbaryon_,\n",
    "                   \"strange_num\":strange_num_,\n",
    "                   \"charge_num\":charge_num_})\n",
    "\n",
    "# 给 pandas 的 DataFrame 数据多添加一列，eta\n",
    "\n",
    "# for fermions (baryons)\n",
    "df.loc[df[\"isbaryon\"]==1, \"eta\"] = +1.0\n",
    "\n",
    "# for bosons (mesons)\n",
    "df.loc[df[\"isbaryon\"]==0, \"eta\"] = -1.0\n",
    "\n",
    "# pdg05.dat 中无反重子，手动添加\n",
    "anti_baryon = df[df['isbaryon']==1].copy()\n",
    "\n",
    "anti_baryon.loc[:, \"pdg\"] = - anti_baryon[\"pdg\"]\n",
    "anti_baryon.loc[:, \"name\"] = \"anti\" + anti_baryon[\"name\"]\n",
    "anti_baryon.loc[:, \"charge_num\"] = -anti_baryon[\"charge_num\"]\n",
    "\n",
    "df = df.append(anti_baryon, ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████| 20/20 [00:20<00:00,  1.29s/it]\n"
     ]
    }
   ],
   "source": [
    "def hrg_eos():\n",
    "    T = np.linspace(0.1, 0.2, 20)\n",
    "    mu = 0  \n",
    "    pressure = []\n",
    "    edensity = []\n",
    "    ndensity = []\n",
    "    cs2 = []\n",
    "    \n",
    "    for Ti in tqdm(T):\n",
    "        pr = 0\n",
    "        ed = 0\n",
    "        nd = 0\n",
    "        for i in range(len(df)):\n",
    "            hadron = df.iloc[i, :]\n",
    "            pr_, ed_, nd_ = eos(Ti, mu, hadron)\n",
    "            pr += pr_\n",
    "            ed += ed_\n",
    "            nd += nd_\n",
    "            \n",
    "        pressure.append(pr)\n",
    "        edensity.append(ed)\n",
    "        ndensity.append(nd)\n",
    "    \n",
    "    cs2 = np.gradient(pressure, edensity)\n",
    "    \n",
    "    return T, np.array(pressure), np.array(edensity), np.array(ndensity), cs2\n",
    "\n",
    "T, Pr, Ed, Nh, Cs2 = hrg_eos()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "论文里的数据点我没有去采集，使用的李文汐同学的数据，特此鸣谢。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAANgElEQVR4nO3ccYjfd33H8efLxE6mtY7lBEmi7Vi6Gsqg7ug6hFnRjbR/JP8USaC4SmnArQ5mETocKvWvKUMQsmm2iVPQWv1DD4nkD1fpECO50lmalMAtOnNE6Fm7/lO0Znvvj99P77hcct/e/e4u3vv5gMDv+/t9fr9758PdM798f/f7paqQJG1/r9rqASRJm8PgS1ITBl+SmjD4ktSEwZekJgy+JDWxavCTfC7Jc0meucLtSfLpJHNJnk7ytsmPKUlaryHP8D8PHLjK7XcB+8Z/jgL/tP6xJEmTtmrwq+oJ4GdXWXII+EKNnALekORNkxpQkjQZOyfwGLuBC0uO58fX/WT5wiRHGf0vgNe+9rV/dMstt0zgy0tSH08++eRPq2pqLfedRPCzwnUrfl5DVR0HjgNMT0/X7OzsBL68JPWR5L/Xet9J/JbOPLB3yfEe4OIEHleSNEGTCP4M8N7xb+vcAbxYVZedzpEkba1VT+kk+TJwJ7AryTzwUeDVAFX1GeAEcDcwB7wEvG+jhpUkrd2qwa+qI6vcXsBfTWwiSdKG8J22ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNTEo+EkOJDmXZC7Jwyvc/uYkjyd5KsnTSe6e/KiSpPVYNfhJdgDHgLuA/cCRJPuXLfs74LGqug04DPzjpAeVJK3PkGf4twNzVXW+ql4GHgUOLVtTwOvHl28ALk5uREnSJAwJ/m7gwpLj+fF1S30MuDfJPHAC+MBKD5TkaJLZJLMLCwtrGFeStFZDgp8Vrqtlx0eAz1fVHuBu4ItJLnvsqjpeVdNVNT01NfXKp5UkrdmQ4M8De5cc7+HyUzb3A48BVNX3gNcAuyYxoCRpMoYE/zSwL8lNSa5j9KLszLI1PwbeBZDkrYyC7zkbSbqGrBr8qroEPAicBJ5l9Ns4Z5I8kuTgeNlDwANJfgB8Gbivqpaf9pEkbaGdQxZV1QlGL8Yuve4jSy6fBd4+2dEkSZPkO20lqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0MCn6SA0nOJZlL8vAV1rwnydkkZ5J8abJjSpLWa+dqC5LsAI4BfwbMA6eTzFTV2SVr9gF/C7y9ql5I8saNGliStDZDnuHfDsxV1fmqehl4FDi0bM0DwLGqegGgqp6b7JiSpPUaEvzdwIUlx/Pj65a6Gbg5yXeTnEpyYKUHSnI0yWyS2YWFhbVNLElakyHBzwrX1bLjncA+4E7gCPAvSd5w2Z2qjlfVdFVNT01NvdJZJUnrMCT488DeJcd7gIsrrPlGVf2yqn4InGP0D4Ak6RoxJPingX1JbkpyHXAYmFm25uvAOwGS7GJ0iuf8JAeVJK3PqsGvqkvAg8BJ4Fngsao6k+SRJAfHy04Czyc5CzwOfKiqnt+ooSVJr1yqlp+O3xzT09M1Ozu7JV9bkn5TJXmyqqbXcl/faStJTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITg4Kf5ECSc0nmkjx8lXX3JKkk05MbUZI0CasGP8kO4BhwF7AfOJJk/wrrrgf+Gvj+pIeUJK3fkGf4twNzVXW+ql4GHgUOrbDu48AngJ9PcD5J0oQMCf5u4MKS4/nxdb+W5DZgb1V982oPlORoktkkswsLC694WEnS2g0Jfla4rn59Y/Iq4FPAQ6s9UFUdr6rpqpqempoaPqUkad2GBH8e2LvkeA9wccnx9cCtwHeS/Ai4A5jxhVtJurYMCf5pYF+Sm5JcBxwGZn51Y1W9WFW7qurGqroROAUcrKrZDZlYkrQmqwa/qi4BDwIngWeBx6rqTJJHkhzc6AElSZOxc8iiqjoBnFh23UeusPbO9Y8lSZo032krSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWpiUPCTHEhyLslckodXuP2DSc4meTrJt5O8ZfKjSpLWY9XgJ9kBHAPuAvYDR5LsX7bsKWC6qv4Q+BrwiUkPKklanyHP8G8H5qrqfFW9DDwKHFq6oKoer6qXxoengD2THVOStF5Dgr8buLDkeH583ZXcD3xrpRuSHE0ym2R2YWFh+JSSpHUbEvyscF2tuDC5F5gGPrnS7VV1vKqmq2p6ampq+JSSpHXbOWDNPLB3yfEe4OLyRUneDXwYeEdV/WIy40mSJmXIM/zTwL4kNyW5DjgMzCxdkOQ24LPAwap6bvJjSpLWa9XgV9Ul4EHgJPAs8FhVnUnySJKD42WfBF4HfDXJfyaZucLDSZK2yJBTOlTVCeDEsus+suTyuyc8lyRpwnynrSQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0MCn6SA0nOJZlL8vAKt/9Wkq+Mb/9+khsnPagkaX1WDX6SHcAx4C5gP3Akyf5ly+4HXqiq3wc+Bfz9pAeVJK3PkGf4twNzVXW+ql4GHgUOLVtzCPi38eWvAe9KksmNKUlar50D1uwGLiw5ngf++EprqupSkheB3wV+unRRkqPA0fHhL5I8s5aht6FdLNurxtyLRe7FIvdi0R+s9Y5Dgr/SM/Vawxqq6jhwHCDJbFVND/j62557sci9WOReLHIvFiWZXet9h5zSmQf2LjneA1y80pokO4EbgJ+tdShJ0uQNCf5pYF+Sm5JcBxwGZpatmQH+Ynz5HuDfq+qyZ/iSpK2z6imd8Tn5B4GTwA7gc1V1JskjwGxVzQD/CnwxyRyjZ/aHB3zt4+uYe7txLxa5F4vci0XuxaI170V8Ii5JPfhOW0lqwuBLUhMbHnw/lmHRgL34YJKzSZ5O8u0kb9mKOTfDanuxZN09SSrJtv2VvCF7keQ94++NM0m+tNkzbpYBPyNvTvJ4kqfGPyd3b8WcGy3J55I8d6X3KmXk0+N9ejrJ2wY9cFVt2B9GL/L+F/B7wHXAD4D9y9b8JfCZ8eXDwFc2cqat+jNwL94J/Pb48vs778V43fXAE8ApYHqr597C74t9wFPA74yP37jVc2/hXhwH3j++vB/40VbPvUF78afA24BnrnD73cC3GL0H6g7g+0Med6Of4fuxDItW3YuqeryqXhofnmL0noftaMj3BcDHgU8AP9/M4TbZkL14ADhWVS8AVNVzmzzjZhmyFwW8fnz5Bi5/T9C2UFVPcPX3Mh0CvlAjp4A3JHnTao+70cFf6WMZdl9pTVVdAn71sQzbzZC9WOp+Rv+Cb0er7kWS24C9VfXNzRxsCwz5vrgZuDnJd5OcSnJg06bbXEP24mPAvUnmgRPABzZntGvOK+0JMOyjFdZjYh/LsA0M/nsmuReYBt6xoRNtnavuRZJXMfrU1fs2a6AtNOT7Yiej0zp3Mvpf338kubWq/meDZ9tsQ/biCPD5qvqHJH/C6P0/t1bV/238eNeUNXVzo5/h+7EMi4bsBUneDXwYOFhVv9ik2TbbantxPXAr8J0kP2J0jnJmm75wO/Rn5BtV9cuq+iFwjtE/ANvNkL24H3gMoKq+B7yG0QerdTOoJ8ttdPD9WIZFq+7F+DTGZxnFfruep4VV9qKqXqyqXVV1Y1XdyOj1jINVteYPjbqGDfkZ+TqjF/RJsovRKZ7zmzrl5hiyFz8G3gWQ5K2Mgr+wqVNeG2aA945/W+cO4MWq+slqd9rQUzq1cR/L8Btn4F58Engd8NXx69Y/rqqDWzb0Bhm4Fy0M3IuTwJ8nOQv8L/Chqnp+66beGAP34iHgn5P8DaNTGPdtxyeISb7M6BTervHrFR8FXg1QVZ9h9PrF3cAc8BLwvkGPuw33SpK0At9pK0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDXx/4aZaro1YsjCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#initialize figure\n",
    "fig, ax = plt.subplots()\n",
    "# fig.set_size_inches(8, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AddLabelAndLegend(ax, xlabel, ylabel):\n",
    "    ax.set_xlabel(xlabel)\n",
    "    ax.set_ylabel(ylabel)\n",
    "    ax.legend()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEGCAYAAACdJRn3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZQU5b3/8feXAUQgJChEkUVGRAFFUZqJS9yIC3oN6A2ukAgxjsrhJh5jFCPGG4y5UaI/43IVjMagKHEPbheXqFGv6AwIyKKCHJYBVES5ShBhnO/vj6fHaZoGumamu2a6P69z+nTXU1U935qlP1P1VD1l7o6IiEgULeIuQEREmh+Fh4iIRKbwEBGRyBQeIiISmcJDREQiaxl3AfnQqVMn79mzZ9xliIg0K7NmzfrE3TtnmlcU4dGzZ08qKyvjLkNEpFkxs+XbmxfbYSszG2Jm75nZEjMbl2H+pWa20MzmmdmLZrZ32vwOZrbKzG7LX9UiIgIxhYeZlQC3AycD/YBzzKxf2mJvAwl3Pwh4BLghbf61wCu5rlVERLYV155HGbDE3Ze6+2ZgGjAsdQF3f8ndNyYnZwLdaueZ2UBgD+C5PNUrIiIp4urz6AqsTJmuAr63g+XPB54FMLMWwI3Aj4Ef1LeALVu2UFVVxaZNm+r7Fs1WmzZt6NatG61atYq7FBFppuIKD8vQlnGQLTMbCSSAY5JNY4Bn3H2lWaa3+Wa9cqAcoEePHtvMr6qq4lvf+hY9e/ZkR+9TaNyddevWUVVVRWlpadzliEgzFVd4VAHdU6a7AavTFzKz44GrgGPc/atk8+HAUWY2BmgPtDazDe6+Vae7u08GJgMkEoltgmnTpk1FFxwAZsbuu+/O2rVr4y5FRJqxuMKjAuhtZqXAKuBs4NzUBczsEGASMMTdP65td/cRKcuMInSqb3O2VjaKLThqFet2i0jjiaXD3N2rgbHADGAR8JC7LzCzCWY2NLnYRMKexcNmNsfMpsdRq4iIbCu2iwTd/RngmbS236S8Pj6L97gXuLexa8uXZcuWceqppzJ//vy4SxERiURjWzVx1dXVcZcgIrINhUfMvv76ay644AIOOOAATjzxRL788kuOPfZYfv3rX3PMMcfwpz/9iQ8++IDDDjuMQYMG8Zvf/Ib27dvHXbaIFLmiGNtqpy65BObMadz3HDAAbr55p4stXryYBx98kLvuuoszzzyTRx99FID169fzyivhAvpTTz2VX/ziF5xzzjnceeedjVuniEg9aM8jZqWlpQwYMACAgQMHsmzZMgDOOuusb5Z54403OOOMMwA499xzt3kPEZF8054HZLWHkCu77LLLN69LSkr48ssvAWjXrl1cJYmI7JT2PJqBww477JvDWdOmTYu5GhERhUezcPPNN3PTTTdRVlbGmjVr+Pa3vx13SSJS5HTYKkY9e/bc6hqPyy67LONyXbt2ZebMmZgZ06ZNI5FI5KtEEZGMFB7NwKxZsxg7dizuzne+8x3uueeeuEsSkSKn8GgGjjrqKObOnRt3GSIi31Cfh4iIRKbwEBGRyBQeIiISmcJDREQiU3jEaHsDHE6ePJk+ffrQp08fEokEL7/88jfztmzZwrhx4+jduzcHHnggZWVlPPvss0A49bd///7079+ffv36MX78eL766quMX0NEpCEUHk3MU089xaRJk3jttdd49913mTx5MiNHjmTVqlUAXH311axZs4b58+czf/58nnzySb744otv1n/ppZd45513eOutt1i6dCnl5eVxbYqIFDCFR5amAj0J37CeyelcuP7665k4cSKdOnUC4NBDD2X06NHcfvvtbNy4kbvuuotbb731mzGx9thjD84888xt3qd9+/bceeedPPHEE3z66ac5qlZEipXCIwtTgXJgOeDJ53JyEyALFixg4MCBW7UlEgkWLlzIkiVL6NGjBx06dMjqvTp06EBpaSmLFy/OQaUi0lD5+qc0FxQeWbgK2JjWtjHZng/uHsu6IpI7+fynNBcUHllYEbG9Ifr168esWbO2aps9ezaJRIJ9992XFStWbNXHsSNffPEFy5YtY7/99stBpSLSEHH/U9pQCo8s9IjY3hCXX345V1xxBevWrQNgzpw5PP7441x44YW0bduW888/n5///Ods3rwZgDVr1nD//fdv8z4bNmxgzJgxnHbaaXTs2DEHlYpIQ+Tzn9Jc0NhWWbiOsDuZ+l9C22R7Q2zcuJFu3bp9M33ppZdy6aWXsnr1ao488kiqq6v58MMPmTt3Lp07dwbgd7/7HePHj6dfv360adOGdu3aMWHChG/e47jjjsPdqamp4fTTT+fqq69uYJUikgs9CIeqMrU3BxbXMXEzGwL8CSgB/uzuf0ibfynwM6AaWAv81N2Xm9kA4A6gA/A1cJ27/21HXyuRSHhlZeVWbYsWLaJv375Z1zuVsDu5gvDDvQ4YkfXa9VNdXc3o0aOpqanh/vvvx8wa7b2jbr+INK7aPo/0f0onk/vPlmyZ2Sx3z3gPiFj2PMysBLgdOAGoAirMbLq7L0xZ7G0g4e4bzexi4AbgLML3+ifuvtjM9gJmmdkMd1+fy5pHkP8faMuWLbnvvvvy/FVFJB9qP0/y/U9pY4mrz6MMWOLuS919MzANGJa6gLu/5O61oTwT6JZsf9/dFydfrwY+BjrnrXIRkUYyAlgG1CSfm0twQHzh0RVYmTJdlWzbnvOBZ9MbzawMaA18kGFeuZlVmlnl2rVrM75psZ7GWqzbLSKNJ67wyHTwPuMnmpmNBBLAxLT2LsB9wGh3r9nmzdwnu3vC3RO1nc2p2rRpw7p164rug9TdWbduHW3atIm7FBFpxuI626oK6J4y3Q1Ynb6QmR1POCR4jLt/ldLeAXgaGO/uM+tTQLdu3aiqqmJ7eyWFrE2bNlud5SUiElVc4VEB9DazUmAVcDZwbuoCZnYIMAkY4u4fp7S3Bh4Hprj7w/UtoFWrVpSWltZ3dRGRohbLYSt3rwbGAjOARcBD7r7AzCaY2dDkYhOB9sDDZjbHzKYn288EjgZGJdvnJE/fFRGRPIntOo98ynSdh4iI7NiOrvPQ8CQiIhKZwkNERCJTeIiISGQKDxERiUzhISIikSk8REQkMoWHiIhEpvAQEZHIFB4iIhKZwkNERCJTeIiISGQKDxERiUzhISIikSk8REQkMoWHiIhEpvAQEZHIFB5xKoIbcYlIYVJ4xGXzZujXD8aOhZkzFSQi0qwoPOKyfj0cfDDcfTccfjjstx/89rfwwQdxVyYislMKj7h897swbRp89BHccw/06BHCY9994Ygj4I47YN26uKsUEclI4RG3Dh1g9Gh48UVYvhyuvx6++ALGjIEuXWDYMHjkEdi0Ke5KRUS+EVt4mNkQM3vPzJaY2bgM8y81s4VmNs/MXjSzvVPmnWdmi5OP8/JbeQ517w6XXw7z5sGcOfCLX0BFBZxxBuy5J/zsZ/DKK1BTE3elIlLkzGPoqDWzEuB94ASgCqgAznH3hSnLHAe86e4bzexi4Fh3P8vMdgMqgQTgwCxgoLt/tr2vl0gkvLKyMncblEtffw0vvQT33QePPgr/+lfoH3nkEejfP+7qRKSAmdksd09kmhfXnkcZsMTdl7r7ZmAaMCx1AXd/yd03JidnAt2Sr08Cnnf3T5OB8TwwJE91519JCRx/PPz1r6F/5P77w2GtI46AJ5+MuzoRKVJxhUdXYGXKdFWybXvOB56Nsq6ZlZtZpZlVrl27toHlNhHt2sGIEeFQ1v77h/6QG27Qab4ikndxhYdlaMv4CWhmIwmHqCZGWdfdJ7t7wt0TnTt3rnehTVLXrvDPf4a+kCuugFGj4Kuv4q5KRIpIXOFRBXRPme4GrE5fyMyOB64Chrr7V1HWLXht24ZTff/zP2HKFBg8GD7+OO6qRKRIxBUeFUBvMys1s9bA2cD01AXM7BBgEiE4Uj8VZwAnmllHM+sInJhsKz5mcM018NBD8PbbMGhQOFNLRCTHYgkPd68GxhI+9BcBD7n7AjObYGZDk4tNBNoDD5vZHDObnlz3U+BaQgBVABOSbcXrjDPg1Vehujp0pP/973FXJCIFLpZTdfOtWZ+qG8Xq1XDaaVBZCb//fegPsUxdRCIiO9cUT9WVXNhrr3AR4VlnwZVXwnnn6cp0EckJhUeh2XVXeOABuPbacGHhccfBhx/GXZWIFBiFRyEyg/Hjw1Xo8+ZBWVkY7kREpJEoPArZj34Er70WLiI88sgwvImISCNQeBS6Qw6Bt96Cgw6C4cPhqqvCeFkiIg2g8CgGXbrAyy/DBReEs7BOPRU+2+44kiIiO6XwiMlUoCfhB9AzOZ1Tu+wCkyfDpEnh3iGDBsH8+bn+qiJSoBQeMZgKlAPLCYNyLU9O5zxAAMrLw+m8GzfC974Xrk4XEYlI4RGDq4CNaW0bk+15cfjhMGsWDBgQrgm54gr1g4hIJAqPGKyI2J4TXbqEm0xdfHEY1v3kk3XPdBHJmsIjBj0itudM69bw3/8Nf/5zOJQ1aBDMnZvvKkSkGVJ4xOA6oG1aW9tkeyzOPz8MrLh5czik9eCDcVUiIs2EwiMGI4DJwN6EO1vtnZweEWdRZWWhHySRgHPPhcsuC6P0iohkoPCIyQhgGVCTfI41OGrtsUc4jfc//gNuvBFOOgk++STuqkSkCVJ4yNZatYJbboF774XXX4eBA+HNN+OuSkSaGIWHZHbeeWFcrBYt4PvfD3siRXDvFxHJjsJDti+RgNmz4Yc/DH0gQ4fqdF4RARQesjMdO4bReG+9FZ57LlxY+PrrcVclIjFTeMjOmcHYsfC//xvGyDrmGPjDH6CmJu7KRCQmCg/J3sCB4TDW8OHhNrennAIffxx3VSISA4WHRNOhQ7iIcNKkMMz7gAHh6nSRiPI+srQ0qtjCw8yGmNl7ZrbEzMZlmH+0mc02s2ozG5427wYzW2Bmi8zsFjOz/FVe2LL6gzYLo/O+9VYIk8GDYcIEDa4oWYt1ZGlpFLGEh5mVALcDJwP9gHPMrF/aYiuAUcADaeseARwJHAQcCAwCjslxyUUh8h/0QQdBZWW4Iv2aa+DEE+HDD/NTrDRrsY8sLQ0W155HGbDE3Ze6+2ZgGjAsdQF3X+bu8wgXYW81C2gDtAZ2AVoBH+W+5MJXrz/o9u1hyhS45x544w04+GB44YWc1SiFoUmMLC0NEld4dAVWpkxXJdt2yt3fAF4C1iQfM9x9UfpyZlZuZpVmVrl27dpGKLnw1fsP2gxGj4aKCujUKeyBXHklfPVVI1cohaLJjCwt9RZXeGTqo8jq8mUz2xfoC3QjBM5gMzt6mzdzn+zuCXdPdO7cuUHFFosG/0EfcEAIkPPPD6fylpVpiHfJqMmNLC2RxRUeVUD3lOluwOos1z0dmOnuG9x9A/AscFgj11eUGuUPum1buOsuePLJcBrvoEFw3XUaoVe20iRHlpZI4gqPCqC3mZWaWWvgbGB6luuuAI4xs5Zm1orQWb7NYSuJrlH/oE89FebPhx/9CMaPhyOOgEX6MUmdJjmytGQtlvBw92pgLDCD8MH/kLsvMLMJZjYUwMwGmVkVcAYwycwWJFd/BPgAeAeYC8x19yfzvhEFqlH/oHffPVwT8tBDsHQpHHII3HSTTukVKQDmRTBSaiKR8MrKyrjLKG4ffggXXgjTp8NRR8Ff/gK9esVdlYjsgJnNcvdEpnm6wlzyY8894Yknwn1C5s4Np/TeeaeGeRdpphQekj9m4T4h8+eHPpCLL4YhQ2Dlyp2vKyJNisJD8q97d5gxA+64Iwzv3r9/uNBQeyEizYbCQ+JhBhddFA5h9e8f9khOOw1WrYq7MhHJwk7Dw8xOMLO7zGxAcro892VJ0ejVK4zOe+ONYW+kT59wRtaWLXFXJiI7kM2exxjgV8BIMxsMDEidmWzDzLokBzwUiaakBC69FBYsgKOPhl/+Mtw75LXX4q5MRLYjm/BY6+7r3f0y4ETCKLaphphZN+BO4P81doFSRHr1gqeegsceg/Xrwym9o0frhlMiTVA24fF07Qt3HwdMSZv/HeAK4HJgU+OVJkXJDE4/PVyNPm4cTJ0K++8fOtd1caFIk7HT8HD3v5vZ1Wb2y+T0rWmLTAD+7u7vse3w6SL1064d/Nd/hQ71Qw6BMWPgsMPC/UNEJHbZnm31Y+CO9EYz+xnwY3d/Ab7ZMxFpPH37wosvwgMPQFVVGKl3zBj47LO4KxMpatmGx5funn6fIID7gJGNWI/ItszgnHPg3Xfh5z8P90/ff3/46191bYhITLIODzPrkt7o7l8BGmtb8uPb34abb4bZs6F3bxg1KpydNW9e3JWJFJ1sw+NG4O9mtndqo5l9F/VzSL4dfDC8+ircfXfoWB8wIATJ8uVxVyZSNLIKD3d/GLgdmGVmT5nZ78zs98DrwB9zWaBIRi1awE9/Cu+/D7/6Ffztb7DffnDJJaDbDovkXNbDk7j7X4FS4CGgFeG03HPcfWqOahPZud12g+uvh8WLwxAnt90G++wDv/0tfPFF3NWJFCzdz0MKy7vvwtVXwyOPQKdO4S6GF10Eu+wSd2UizY7u5yHFo08fePhheOut0DdyySXhzKwpU3SRoUgjUnhIYRo0CF54AZ5/PuyBnHdeCJPp03V6r0gjUHhIYTv+eKioCHsjW7bAsGHw/e+Hs7VEpN4UHlL4zGD48DBq7+TJsGxZuD7kuOPgmWe0JyJSDwoPKR4tW8IFF8CSJeH+IUuWwL/9Gxx0ULhaffPmuCsUaTYUHlJ8dt013D9k6dLQkW4WLjLcZx+YOBH+7//irrDJmAr0JHxQ9ExOi0CM4WFmQ8zsPTNbYmbbDKhoZkeb2Wwzqzaz4WnzepjZc2a2yMwWmlnPfNUtBaRVK/jxj8PIvf/zP+GsrMsvhx49wnOR3xJ3KlAOLAc8+VyOAkSCWMIjecfB24GTgX7AOWbWL22xFcAo4IEMbzEFmOjufYEyQHcLkvozg5NOCqP3VlbCKaeEw1qlpWGPZP78uCuMxVVA+mioG5PtInHteZQBS9x9qbtvBqYBw1IXcPdl7j6PtLGzkiHT0t2fTy63YTsj/opEN3AgPPhg6A+56KJwllb//qFv5OWXi6pzfUXEdikucYVHV2BlynRVsi0b+wHrzewxM3vbzCZmune6mZWbWaWZVa7VWEcSVWkp3HILrFgB114bTvc97rhw/cikSfD553FXmHM9IrZLcYkrPCxDW7b/0rUEjgIuI9xPfR/C4a2t38x9srsn3D3RuXPn+tYpxW733cMQJ8uXw513hjOyLroI9twzHNJ69dWC3Ru5Dmib1tY22S4SV3hUAd1TprsBqyOs+3bykFc18ARwaCPXJ7K1XXeFCy8MnetvvRU62h97LFwv0qdPGJzxww/jrrJRjQAmA3sT/tvbOzk9Is6ipMmIKzwqgN5mVmpmrYGzgekR1u1oZrW7E4OBhTmoUWRbZnWHrtasgXvvhT32gHHjoFu3cAX7k09CdWHcI20EsIzQ8bgMBYfUiSU8knsMY4EZwCLgIXdfYGYTzGwogJkNMrMq4AxgkpktSK77NeGQ1Ytm9g7hn6K74tgOKXLt2oUxs/75zzCa7y9/CW++CUOHhtN9r7wyDBUvUoA0JLtIY9qyJQx5cvfd8PTTUFMTDm2NHBn2Sr773bgrFMmahmQXyZdWrUJITJ8OK1fC738fDm+Vl0OXLnDsseEsrpUrd/pWIk2Z9jxEcs0d5s0LHeyPPVZ30WFZGfz7v4dH797x1iiSwY72PBQeIvn2/vt1QVJREdoOPLAuSA46KHTMi8RM4aHwkKZqxQp4/PEQJLXXjPTqFULk9NPDmV0tW8ZdpRQphYfCQ5qDjz4KfSWPPhrG2aquhg4dQj/J4MHwgx/AAQdor0TyRuGh8JDmZv16mDEjhMiLL4bh4yFcU1IbJIMHh2FURHJE4aHwkOZu2TL4xz9CkPzjH3VXs5eWhiD5wQ/C2Ft77BFrmVJYFB4KDykk7rBoUd1eycsv193A6sAD4aijQl/JoEHQty+UbDNuqEhWFB4KDylkX38Ns2fXhcmbb8IXX4R57drBoYeGIEkkwnOvXuo3kawoPBQeUkxqasLpwBUVdY85c2DTpjC/Y8e6IKkNla5dFSiyDYWHwkOK3ZYtsGDB1oHyzjthrwXCEPMHHBBGCO7bNzz36QN77aVQKWIKD4WHyLa+/DLskVRWwqxZoR9l0aK6Q14A3/rW1oFS+9yrVxiKRQqawkPhIZId9zAW17vvhiBJfV61qm65li1h333DsCrdu4dHjx51z3vtpXApADsKD126KiJ1zMIH/157hetIUn3+Obz3Xl2gLFoUrj959dVwXUr6+3TpUhcoqeHSvXs4pbhTJ2ibfq9CaS4UHiKSnQ4d6jrZ023YEEYKXrFi6+eVK8OhsSefrOuwT9WmTQiR1Mfuu2du2223cPZYu3ZhPfXFxErhISIN17596A/p2zfzfHf45JO6QFm7NkyvWxeea18vXx5ef/bZjr9eixZ1QdK+fd3rTNO77BIOs0V5lJSEkwlqasKj9nX6c3pblEd1dd3rK69sdiMrKzxEJPfMoHPn8Dj00J0vX10dAiQ1WD79FP71r/DYsCHz688/D302qe1btoRHXEpKtn3UBlTt48IL46uvnhQeItL0tGxZFzaNpaYmhFI2j6+/Dns3JSXhOfV1+nN6W+qjReHeb0/hISLFoUULaN06PKTBCjcWRUQkZxQeIiISWWzhYWZDzOw9M1tiZuMyzD/azGabWbWZDc8wv4OZrTKz2/JTscRhKtCT8IvaMzktmel7JfkUS3iYWQlwO3Ay0A84x8z6pS22AhgFPLCdt7kWeCVXNUr8pgLlwHLAk8/l6EMxE32vJN/i2vMoA5a4+1J33wxMA4alLuDuy9x9HlCTvrKZDQT2AJ7LR7ESj6uAjWltG5PtsjV9ryTf4gqPrsDKlOmqZNtOmVkL4EbgVztZrtzMKs2scu3atfUuVOKzImJ7MdP3SvItrvDINK5AtiM0jgGecfeVO1rI3Se7e8LdE50b81xxyZseEduLmb5Xkm9xhUcV0D1luhuwOst1DwfGmtky4I/AT8zsD41bnjQF1wHpw+a1TbbL1vS9knyL6yLBCqC3mZUCq4CzgXOzWdHdR9S+NrNRQMLdtzlbS5q/2h/0VYTDLz0IH4YjtrtG8dL3SvItlvBw92ozGwvMAEqAe9x9gZlNACrdfbqZDQIeBzoCPzSz37r7AXHUK/EZgT4As6XvleSTbgYlIiIZ7ehmULrCXEREIlN4iIhIZAoPERGJTOEhIiKRKTxERCQyhYeIiESm8BBBw5mLRKXb0ErRqx3OvHZU2trhzEEX3Ylsj/Y8pOhpOHOR6BQeUvQ0nLlIdAoPKXoazlwkOoWHFD0NZy4SncJDit4IYDKwN+EuZXsnp9VZLrJ9OttKBA1nLhKV9jxERCQyhYeIiESm8BARkcgUHiIiEpnCQ0REIlN4iIhIZAoPERGJLLbwMLMhZvaemS0xs3EZ5h9tZrPNrNrMhqe0DzCzN8xsgZnNM7Oz8lu5SHbqM8y7hoaX5iKWiwTNrAS4HTgBqAIqzGy6uy9MWWwFMAq4LG31jcBP3H2xme0FzDKzGe6+Pg+li2SlPsO8a2h4aU7i2vMoA5a4+1J33wxMA4alLuDuy9x9HlCT1v6+uy9Ovl4NfAx0zk/ZItmpzzDvGhpempO4wqMrsDJluirZFomZlQGtgQ8yzCs3s0ozq1y7dm29CxWpj/oM866h4aU5iSs8LEObR3oDsy7AfcBod69Jn+/uk9094e6Jzp21YyL5VZ9h3jU0vDQncYVHFdA9ZbobsDrblc2sA/A0MN7dZzZybSINVp9h3jU0vDQncYVHBdDbzErNrDVwNjA9mxWTyz8OTHH3h3NYo0i91WeYdw0NL82JuUc6WtR4X9jsFOBmoAS4x92vM7MJQKW7TzezQYSQ6AhsAj509wPMbCTwF2BBytuNcvc52/taiUTCKysrc7YtIiKFyMxmuXsi47y4wiOfFB4iItHtKDx0hbmIiESm8BARkcgUHiIiEpnCQ0REIlN4iIhIZAoPERGJTOEhIiKRKTxERCQyhYeIiESm8BARkcgUHiIiEpnCQ0REIlN4iIhIZAoPERGJTOEhIiKRKTxERCQyhYeIiESm8BARkcgUHiIiEpnCQ0REIlN4iIhIZLGFh5kNMbP3zGyJmY3LMP9oM5ttZtVmNjxt3nlmtjj5OC9/VYuICMQUHmZWAtwOnAz0A84xs35pi60ARgEPpK27G3AN8D2gDLjGzDrmumYREakT155HGbDE3Ze6+2ZgGjAsdQF3X+bu84CatHVPAp5390/d/TPgeWBIPooWEZEgrvDoCqxMma5KtuV6XRERaQQtY/q6lqHNG3NdMysHypOTG8zsvSzfP986AZ/EXUQjKaRtgcLankLaFiis7WnK27L39mbEFR5VQPeU6W7A6gjrHpu27svpC7n7ZGBy/crLHzOrdPdE3HU0hkLaFiis7SmkbYHC2p7mui1xHbaqAHqbWamZtQbOBqZnue4M4EQz65jsKD8x2SYiInkSS3i4ezUwlvChvwh4yN0XmNkEMxsKYGaDzKwKOAOYZGYLkut+ClxLCKAKYEKyTURE8iSuw1a4+zPAM2ltv0l5XUE4JJVp3XuAe3JaYP40+UNrERTStkBhbU8hbQsU1vY0y20x92z7qUVERAINTyIiIpEpPEREJDKFR44U2thd9d0eMxtgZm+Y2QIzm2dmZ+W38m015GeTnN/BzFaZ2W35qXjHGvi71sPMnjOzRWa20Mx65qvuTBq4LTckf88WmdktZpbpmrC8ymJ7Lk1+3+eZ2YtmtnfKvCb3ObAVd9ejkR9ACfABsA/QGpgL9EtbpidwEDAFGJ7SvhuwNPncMfm6YzPenv2A3snXewFrgO80x21Jmf8nwphrtzXn37XkvJeBE5Kv2wNtm+O2AEcAryffowR4Azi2Gfxsjqv9ngMXA39Lvm5ynwPpD+155Eahjd1V7+1x9/fdfXHy9WrgY6BzfsrOqCE/G8xsILAH8Fw+is1CvbcnORhpS3d/PrncBnffmKe6M2nIz/QwCyMAAAKkSURBVMaBNoQP6V2AVsBHuS95h7LZnpdSvuczqTvDtCl+DmxF4ZEbhTZ2V6PUZGZlhD/uDxqprvqo97aYWQvgRuBXOairvhrys9kPWG9mj5nZ22Y2MTnidVzqvS3u/gbwEmHPdg0ww90XNXqF0UTdnvOBZ+u5bt4pPHIj52N35VmDazKzLsB9wGh33+Y/+jxqyLaMAZ5x95U7XTJ/GrI9LYGjgMuAQYTDK6Map6x6qfe2mNm+QF/Cf+5dgcFmdnQj1lYfWW+PmY0EEsDEqOvGReGRGw0du6u+6+ZKg2oysw7A08B4d5/ZyLVF1ZBtORwYa2bLgD8CPzGzPzRueZE19Hft7eRhlWrgCeDQRq4vioZsy+nAzOShtw2E/+APa+T6ospqe8zseOAqYKi7fxVl3TgpPHKj0Mbuqvf2JJd/HJji7g/nsMZs1Xtb3H2Eu/dw956E/9anuPs2Z9DkWUN+1yqAjmZW2wc1GFiYgxqz1ZBtWQEcY2YtzawVcAxh6KM47XR7zOwQYBIhOD5OmdUUPwe2FnePfaE+gFOA9wnH969Ktk0g/JJAOExQBfwLWAcsSFn3p8CS5GN03NvSkO0BRgJbgDkpjwHNcVvS3mMUTeBsq0b4XTsBmAe8A9wLtG6O20I4s2kSITAWAjfF/XPJcnteIHTs1/5tTE9Zt8l9DqQ+NDyJiIhEpsNWIiISmcJDREQiU3iIiEhkCg8REYlM4SEiIpHFdidBkWJmZrsDLyYn9wS+BtYmp8s8jIUk0mTpVF2RmJnZfwIb3P2Pcdciki0dthIRkcgUHiIiEpnCQ0REIlN4iIhIZAoPERGJTOEhIiKR6VRdERGJTHseIiISmcJDREQiU3iIiEhkCg8REYlM4SEiIpEpPEREJDKFh4iIRPb/ASTJaulI3LdgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Cs2\n",
    "ax.cla()\n",
    "data = pd.read_csv('./data/C^2_s_T.csv')\n",
    "T_lqcd = data['T'] / 1000\n",
    "Cs2_lqcd = data['C^2_s']\n",
    "ax.scatter(T_lqcd, Cs2_lqcd, color='cyan', label='LQCD')\n",
    "ax.plot(T, Cs2, color='red', label='hrg')\n",
    "AddLabelAndLegend(ax, 'T', r'$C^2_s$')\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXiU5dn+8e/FLuDGUlwQQkVFRIsSEbUoyitVcUFFikBbEcWWuhUUERBcQAtYFRWF4KvoTxYXikstlB0VcAkqymIFARGJJW+sCrKG3L8/7skYYgJZZuZ5Zub8HEeOmTyznVkmV+7n3sw5h4iICECVoAOIiEh4qCiIiEiUioKIiESpKIiISJSKgoiIRFULOkBFNWjQwGVkZAQdQ0QkaSxbtuz/nHMN93efpC0KGRkZZGdnBx1DRCRpmNmXB7qPTh+JiEiUioKIiESpKIiISFTS9imUZM+ePWzatImdO3cGHSXhatWqRePGjalevXrQUUQkiaVUUdi0aRMHH3wwGRkZmFnQcRLGOUdeXh6bNm2iWbNmQccRkSSWUqePdu7cSf369dOqIACYGfXr10/LFpKIxFZKFQUg7QpCoXT9ukUktlKuKIhIEnAO7rgDVq4MOokUo6IQYxs2bKBVq1ZBxxAJtxkz4KGHQBNQQ0dFIQD5+flBRxAJzt69MGwYnHAC9OwZdBopRkUhDvbu3csNN9zASSedRKdOndixYwcdOnRg8ODBnHvuuYwdO5YvvviCdu3acfrppzNs2DDq1q0bdGyRxJg2zZ82uu8+qJZSAyBTQur+RG67DT7+OLbP2bo1PProAe+2Zs0apk6dysSJE+nWrRvTp08H4LvvvmPRokUAXHLJJdx6661cc801jB8/PrY5RcJqzx645x741a+ga9eg00gJ1FKIg2bNmtG6dWsA2rRpw4YNGwD47W9/G73P0qVLufrqqwHo0aNHwjOKBOK552DtWrj/fqiiPz9hlLothTL8Rx8vNWvWjF6vWrUqO3bsAKBOnTpBRRIJ3q5d/pTRGWfAJZcEnUZKoVIdkHbt2kVPK02bNi3gNCIJkJUFX30FI0aA5tWElopCQB599FEefvhh2rZtS05ODoceemjQkUTi58cfYeRI6NABOnYMOo3sR+qePgpIRkYGK1asiH5+++23l3i/o48+mnfffRczY9q0aWRmZiYqokjijRsH//kPTJ+uVkLIqSgEZNmyZdx000045zjssMN45plngo4kEh/ffw+jRsFFF8HZZwedRg5ARSEg7du3Z/ny5UHHEIm/Rx6Bb7/1I44k9NSnICLxk5cHDz8MV14JbdoEnUbKQEVBROJn9GjYts0PRZWkoKIgIvHxzTfw+OPQowecdFLQaaSMVBREJD4eeAB27/bLWkjSUFGIsdIWtsvKyqJFixa0aNGCzMxMFi5cGL1tz549DBo0iOOOO45WrVrRtm1bZs6cCfghrieffDInn3wyLVu2ZOjQoezatSsRX4pIxW3cCBMmQO/e0Lx50GmkHFQUEuAf//gHEyZM4J133uGzzz4jKyuLXr168fXXXwNw9913k5OTw4oVK1ixYgVvvPEGW7dujT5+wYIFfPrpp7z//vusW7eOvn37BvWliJRN4Uiju+8ONoeUW0KLgpk9Y2ZbzGxFkWP1zGyOma2JXB6eqDyTgQz8NyEj8nk8jBo1ijFjxtCgQQMATjvtNHr37s24cePYvn07EydO5PHHH4+umdSoUSO6dev2s+epW7cu48eP59VXX+Xbb7+NU1qRSlqzBp59Fm68EZo0CTqNlFOiWwqTgAuLHRsEzHPOHQfMi3wed5OBvsCXgItc9iU+hWHlypW0KTYcLzMzk1WrVrF27VqaNGnCIYccUqbnOuSQQ2jWrBlr1qyJQ1KRGLj3XqhRAwYPDjqJVEBCi4Jz7i2g+L+4lwPPRa4/B3RJRJYhwPZix7ZHjieCcy6Qx4rE1YoVMGUK3HwzHHFE0GmkAsLQp9DIOZcDELn8RWl3NLO+ZpZtZtm5ubmVetGN5TxeGS1btmTZsmX7HPvwww/JzMykefPmbNy4cZ8+hP3ZunUrGzZs4Pjjj49DUpFKGj4c6taFgQODTiIVFIaiUGbOuSznXKZzLrNhw4aVeq7SznTG4wzowIEDufPOO8nLywPg448/ZsaMGdx4443Url2bPn36cMstt7B7924AcnJyeOGFF372PNu2baNfv3506dKFww9PWNeLSNm8+y78/e/Qvz/Urx90mtBIVN9lrIRh7aP/mNmRzrkcMzsS2JKIFx2J70MoegqpduR4ZWzfvp3GjRtHP+/fvz/9+/dn8+bNnH322eTn5/PNN9+wfPlyCgvbiBEjGDp0KC1btqRWrVrUqVOH+4rMAD3vvPNwzlFQUMAVV1zB3RrRIWHjnC8GjRpBKSsDp6PCvsvCvzOFfZcAPQNJdGCW6PPTZpYB/MM51yry+Rggzzn3VzMbBNRzzh2w7ZmZmemys7P3ObZ69WpOPPHEMmeZjO9D2IhvIYwk/j+o/Px8evfuTUFBAS+88AIWw2WEy/v1i8TMyy9Dt25+I50bbgg6TWhk4AtBcU2BDQlN4pnZMufcftfpT2hRMLOpQAegAfAfYDjwKvAS/u/yRuBq59wBx1vGoiikmnT/+iUgu3ZBy5ZQuzZ8/DFUrRp0otCogh/dWJwBBQnOAmUrCgk9feScu6aUm7QVk0iyGjcO1q2DWbNUEIppQskthTDP3kiqjuaySNfhmun6dUvA8vL87OXf/MZ/yD5G4vsqi4pF32U8pVRRqFWrFnl5eWn3B9I5R15eHrVq1Qo6iqSb+++HH36Ahx4KOkko9QSy8H0IFrnMIrydzBCO0Ucx07hxYzZt2kRl5zAko1q1au0z6kkk7tas8aeO+vSBVq2CThNaPQl3ESgupYpC9erVadasWdAxRNLDoEFQs6Y20EkxKXX6SEQS5O23/US1O+/UchYpRkVBRMqnoAAGDICjj/aXklJS6vSRiCTAtGnwwQcwaZKfmyApRS0FESm7nTvhrrvg1FPhd78LOo3EgVoKIlJ2Y8f6rTYnTYIq+p8yFemnKiJlk5sLDzwAl14K550XdBqJExUFESmbe+6BH3+E0aODTiJxpKIgIgf22WcwYYLfd7lFi6DTSBypKIjIgQ0cCHXq+NaCpDR1NIvI/i1YAG+8AX/9K1Ryx0MJP7UURKR0e/f6CWpNm8KttwadRhJARUFESpeVBR99BKNGQQxX4U22fYvTiU4fiUjJtmyBwYPh/PP9Vpsxkoz7FqcTtRREpGR33umHoI4bBzHcS3wIPxWEQtsjxyV4Kgoi8nOLF/tZywMGxHwI6sZyHpfEUlEQkX3l50O/fnDMMTB0aMyfvrT9icO8b3E6UVEQkX2NGweffAKPPurnJsRYMu5bnE5UFETkJzk5cPfdcOGFcMUVcXmJZNy3OJ1o9JGI/OSOO2DXLnjssZh2LheXbPsWpxO1FETEW7gQJk/2o46OOy7oNBIQFQURgT174M9/howMv4mOpC2dPhIR36m8ahW8/jocdFDQaSRAaimIpLtNm+Dee/3mOZdeGnQaCZiKgki669/fL3w3dmzQSSQEVBRE0tns2fDyyzBkCDRrFnQaCQEVBZF0tWsX3HQTNG8Ot98edBoJidB0NJvZX4DrAQd8CvR2zu0MNpVICvvb32DNGpg1K6bLYktyC0VLwcyOBm4BMp1zrYCqQPdgU4mksA0bYMQIuOoq+M1vgk4jIRKKohBRDTjIzKrhl0LZHHAekdTkHNx2m5+x/MgjQaeRkAlFUXDOfQ08hF89Nwf43jk3u/j9zKyvmWWbWXZubm6iY4qkhldegddeg+HD/UqoIkWEoiiY2eHA5UAz4Cigjpn1Kn4/51yWcy7TOZfZUBuIi5Rfbq6fuZyZ6YeiihQTiqIA/A+w3jmX65zbA/wdOCvgTCKp59Zb4bvv4JlnoFpoxplIiISlKGwE2plZbTMzoCOwOuBMIqnltddg6lS/cc7JJwedRkIqFEXBOfce8ArwIX44ahX8EusiEgv//S/88Y/wq19pwTvZr9C0H51zw4HhQecQSUn9+/v+hDffhOrVg04jIRaKloKIxNHMmTBpkt8n4bTTgk4jIaeiIJLKfvgB+vaFE0+EYcOCTiNJIDSnj0QkDgYOhM2bYckSqFkz6DSSBNRSEElV8+fDhAnwl7/AGWcEnUaShIqCSCratg2uv97vtXz//UGnkSSi00ciqWjIEFi/Ht56S9trSrmopSCSat55Bx5/3O+V0L590GkkyagoiKSSHTvguuugaVN48MGg00gS0ukjkVQybJjfOGfuXKhbN+g0koTUUhBJFe+9Bw8/DDfcAB07Bp1GkpSKgkgq2L4deveGo46CMWOCTiNJTKePRFLBgAGwejXMng2HHhp0GkliaimIJLtXX4Xx4+H22+GCC4JOI0lORUEkmX39NfTp4xe6Gzky6DSSAlQURJLV3r3wu9/Bzp0wZQrUqBF0IkkB6lMQSVZjxsCCBfD003DCCUGnkRShloJIEpgMZODfsBnArPffh7vvhquv9pPVRGJELQWRkJsM9AW2Rz7P27qV43r04Mcjj6TOhAlgFmA6STVqKYiE3BB+KggAj998Mxnr13Pt5Mlw+OFBxZIUpaIgEnIbi1zvPnUq1z73HCOHDGG6FruTOFBREAm5JpHLjPXrGf/HP7LkzDO5b9iw6HGRWFJREAm5kcDB+flM7tkTgJ6TJ1OzWjU0K0HiQR3NIiHXEzj5/vs5ZelSekyZgmvWjKzIcZFYU1EQCbu33+aUESPg979nyjXXBJ1GUpxOH4mE2X//Cz17QrNm8MQTQaeRNKCWgkhYOef3RsjJgcWL4eCDg04kaUBFQSSsHnoIpk+H0aOhbdug00iaOGBRMLOLC68C1wMTnXP/jGsqkXQ3dy4MGgRdu/olsUUSpCx9CvcBJwENgNqRSxGJlw0boHt3aNECnnlGy1hIQpWlKJwD1AV2Aiudc8/HI4iZHWZmr5jZZ2a22szOjMfriITajh1w1VWwZw/MmKF+BEm4A54+cs5tB4abWUf2XYIl1sYCs5xzXc2sBr5VIpI+nIM//Qk+/BDeeAOOPz7oRJKGytTRbGYPOeduB+bFI4SZHYJvkVwL4JzbDeyOx2uJhNaTT8Jzz8Hw4XDJJUGnkTRV1nkK58c1BfwSyAWeNbOPzOxpM6sT59cUCY/Fi+G226BzZxg2LOg0ksbCMnmtGnAa8JRz7lTgR2BQ8TuZWV8zyzaz7Nzc3ERnFImPzZv9KKOMDHjhBagSlrelpKOy/vb9yszWm9nrZvaAmV1jZiebWfUY5dgEbHLOvRf5/BV8kdiHcy7LOZfpnMts2LBhjF5aJD6K75Y2uaQ77d7td0/butV3LB92WOICipSgrEXhE+Bs4AkgD+gEPAv8n5mtqGwI59w3wFdmVrjRbEdgVWWfVyQohbulfQm4yGVfSigMt90GS5b4oaetWiU2pEgJyjyj2Tm3GdgMzC48ZmYGNI9RlpuByZGRR+uA3jF6XpGEK75bGpHPh1BkddNnn4WnnoI77oBu3RIZT6RUZS0KM8zMnHOu6MHI52tiEcQ59zGQGYvnEgnaxgMdz872w087doQHHkhQKpEDK+vpo6OAZWY2zcyuNbMj4hlKJNmVtitaE4DcXLjySmjUCKZNg2pagkzCo0y/jc65PwKYWQvgImCSmR0KLABmAYudc3vjllIkyYzE9yEUPYVUG3hw1y7fsbxlix+G2kCrxki4lGvsm3PuM+fcI865C/FzF94Brgbe2/8jRdJLTyALaIpfSbIpkFVQwDXXXQeLFvmO5TZtAs0oUpJytVvNbKlz7kwA59wOM3sbyHHO3RyXdCJJrCfFtswcOhSmTPF9CD16BJRKZP/KO0umJoCZPQzgnNsKPBnrUCIpZ8IEePBB6NvXL4ktElLlLQpmZr8AekWGowIcFONMIqnlzTehXz+4+GIYN05LYUuolXfYw134foQpwCNm9jnhWSpDJHyys/0chNat4cUXNdJIQq9cv6HOuVnA8QCR/Q6uBvrEIZdI8tuwwa922rChby3UrRt0IpEDqvC/Lc65pcDSGGYRSR3ffgsXXQS7dsGCBXCEpvZIclBbViTWdu2CK66Adetgzhw48cSgE4mUmYqCSCwVFMAf/gBvvQVTp8I55wSdSKRc1EksEkt33eU7lEeNgu7dg04jUm4qCiKx8uSTMHq0X+jujjuCTiNSISoKIrHw+utw881w6aXw2GOaiyBJS0VBpBRl2jkNYPZsv8hdmza+H0FzESSJ6bdXpASFO6cVrnJauHMaFFvPaP58uPxyP8Jo1iyoUyeBKUViTy0FkRLsb+e0qLff9qeLjj0W5s6FevUSlk8kXlQUREpwwJ3Tli71axk1aQLz5mlfBEkZKgoiJdjvzmkffAAXXuhnKc+b53dQE0kRKgoiJRiJ3ymtqNrAkx99BJ06Qf36vj/hqKMCSCcSPyoKIiUoaee0lz79lIsvuAAOPtgXhGOOCTSjSDxo9JFIKfbZOW3VKujYEWrW9AvcZWQEF0wkjtRSEDmQzz/3BaFKFd9COPbYoBOJxI1aCiL788UXcP75sHcvLFwIJ5wQdCKRuFJRECnNl1/6grBjhz9l1LJl0IlE4k6nj0RKsno1/PrX8MMPfk+EU04JOpFIQqgoiBT33nu+IOzZ41sIp50WdCKRhFFRECnqX//yp4wOOwwWL4bWrYNOJJJQKgoihaZO9WsZHXecLwgaZSRpSEVBBOCJJ6BnTzjzTFi0yC9hIZKGQlUUzKyqmX1kZv8IOoukCedg+HC/Qc5ll/nlrw89NOhUIoEJVVEAbgVWBx1CUk+JG+bs3Qv9+sF998F118Err8BBBwUXUiQEQlMUzKwx0Bl4OugskloKN8z5EnCRy5t27eLL7t1h/HgYNAieflo7pokQrslrjwIDgYNLu4OZ9SWyAVaTJqUtbiyyr+Ib5tTdupVXunSh6fz58Le/Qf/+QUUTCZ1QtBTM7BJgi3Nu2f7u55zLcs5lOucyGzZsmKB0kuyKbpjTcMsWFpx3HucuWsTvn39eBUGkmLC0FM4GLjOzi4FawCFm9oJzrlfAuSQFNMGfMjpl+XJeu/xyfrFlC5e/9horO3cOOppI6ISipeCcu8s519g5lwF0B+arIEisjAR6vvwyS846i2r5+Zy7aBELO3dmZNDBREIoFEVBJG4KCug5dCgvdOvG6tatOT07m9zTTyeLInsliEhUWE4fRTnnFgILA44hqeCHH6BXL3jjDejTh8xx48ipWTPoVCKhFrqiIBITn38OXbr4yyee8PMRzIJOJRJ6KgqSembNgu7d/byDuXOhQ4egE4kkDfUpSOpwDsaMgc6doWlTyM5WQRApJ7UUJDXs2AHXXw9TpkDXrjBpEtSpE3QqkaSjloIkv40boX17XxBGjICXXlJBEKkgtRQkub38MvTt6xe3e+01v9KpiFSYWgqSnLZtgz59oFs3vynOhx+qIIjEgIqCJJ/sbL9v8rPPsmLwYJovXkyV5s1/WhJbRCpMRUGSR0EBjB7td0fbsYM58+dzxsiRfFG9enRJ7L6oMIhUhoqCJIevv4YLLoA774TLL4fly7mhQ4d9lsQGv0T2kCDyiaQIFQUJvxkz4JRT4N13/WY4L78M9ertsyR2UaUdF5EDU1GQ8PrxR7jxRrjySsjI8J3JffpEl6sobZslbb8kUnEqChJOH3wAmZkwcSIMHAhLl8IJJ+xzl5FA7WIPqx05LiIVo6Ig4bJ1K9xyC5xxhl/ldM4cGDUKatT42V17AllAU8Ail1oSW6RyNHlNwuPVV+Gmm2DzZr+q6ciRcOih+31IT1QERGJJLQUJ3qZNcMUV/qNePViyxC93fYCCICKxp6Igwdm7Fx5/HFq2hH/9C/76V1i2DNq1CzqZSNrS6SMJxvLlfs2i99+HTp3gqafgl78MOpVI2lNLQRJr+3Y/Aa1NG1i/HiZP9pviqCCIhIJaCpIYzvmO5AEDfDHo08cvWVGvXtDJRKQItRQk/pYsgV//2k9CO+ggWLTIz0xWQRAJHRUFiZ81a/wuaGefDevWQVaW70s455ygk4lIKVQUJPZyc+Hmm/2oolmz4N57Ye1auOEGqPbTGcvJQAb+lzADrW4qEgbqU5DY2b4dHn3UDy3dvt0XgXvugUaNfnbXyfhlrgtXOS1c9ho0GU0kSGopSOXt3QvPPON3QBsyBM4/H1as8MNMSygI4Je31rLXIuGjoiAVV1AA06fDqaf60UTHHANvveVHGbVosd+HatlrkXBSUZDy27MHJk3yfQZdu8LOnfDSS34l0/bty/QUWvZaJJxUFKTsduzwaxI1bw69e0OtWvDii7B6NVx9dXSfg7LQstci4aSOZjmw77/3/QOPPAJbtvghpk89BRddVK5CUFRhZ/IQ/CmjJviCoE5mkWCpKEjpcnNh7FjfOvj+e7jwQhg8uMyniA5Ey16LhE8oioKZHQM8DxwBFABZzrmxwaZKY2vX+tVLJ070/QVXXQV33QWnnRZ0MhGJs1AUBSAfGOCc+9DMDgaWmdkc59yqoIOljfx8ePNNePJJmD3bTzLr1csvXneAkUQikjpCURScczlATuT6VjNbDRwNqCjEW06OX4coK8tvdnP00XDffXD99XDkkUGnE5EEC0VRKMrMMoBTgfdKuK0vkYmvTZpo8GKFOQcLF/rO4hkzfCuhUyd/yuiSS/ZZikJE0kuohqSaWV1gOnCbc+6H4rc757Kcc5nOucyGDRsmPmCy++47eOwxOOkkP+t43jy49Vb4/HO/81mXLqUWBK1TJJIeQvMvoZlVxxeEyc65vwedJ2Xs2eP/4E+Z4mca79gBZ5wBzz3n5xYcdNABn0LrFImkj1AUBTMz4H+B1c65h4POk/QKCvweBlOm+JnGeXlQvz5ce63vKyjnKKL9rVOkoiCSWkJRFICzgd8Bn5rZx5Fjg51z/wwwU/JZudJvbzllCnz5pW8FdOkCPXv6PoPq1Sv0tFqnSCR9hKIoOOfeASo2NTbdffUVTJ3qi8Enn0DVqr4AjBjhC0LdupV+iSb4U0YlHReR1BKKoiDl4Jz/4//66/4jO9sfb9fOjx7q1g1+8YuYvuRI9u1TAK1TJJKqVBSSwe7dfl/jwkKwcaNfc6hdO3jwQd9hfOyxcXt5rVMkkj5UFMLq229h5kxfBGbOhK1bfR9Bp04wfDh07lzqBjbxoHWKRNKDikJY5OfDBx/4uQNz58I77/gdzY44Arp3h8sug44dyzSEVESkolQUguKc37Jy3jz/sWiRbw0AtG4Ngwb5QpCZCVX2nWM4GZ3KEZH4UFFIpPXrfyoC8+f7vQnAb1rTo4dvCZx3HjRoUOpTaCKZiMSTikK85Of7UUJLl/qJZEuWwIYN/rYjjoALLvBFoGNHKMc6TppIJiLxpKIQK3l5vgAUFoH334ftkT/fRx4JZ50F/fv7InDiiRXesUwTyUQknlQUKmrXLnj+eV8Ali6Ff//bH69a1fcJ9OnjC8GZZ/qWQAWLQHGaSCYi8aSiUFHVqsGAAVCzpv/Df+21/jIzE+rUidvLaiKZiMSTikJFVa0Kn33mTw3FqBVQFppIJiLxpKJQGUcdFcjLaiKZiMRLqDbZERGRYKkoiIhIlIpCBWhrShFJVepTKCfNKBaRVKaWQjntb0axiEiyU1EoJ80oFpFUpqJQTqXNHNaMYhFJBSoK5TQSP4O4KM0oFpFUoaJQTj2BLKApYJHLLNTJLCKpQaOPKkAzikUkVamlICIiUSoKIiISpaIgIiJRKgoiIhKloiAiIlHmnAs6Q4WYWS4l70yZaA2A/ws6xAEkQ0ZQzlhKhoyQHDmTISOULWdT51zD/d0haYtCWJhZtnMuM+gc+5MMGUE5YykZMkJy5EyGjBC7nDp9JCIiUSoKIiISpaJQeVlBByiDZMgIyhlLyZARkiNnMmSEGOVUn4KIiESppSAiIlEqCiIiEqWiUAozu9DM/m1ma81sUAm3n2NmH5pZvpl1LXbbH8xsTeTjD2HMaWatzWypma00s0/M7Ldhy1jk9kPM7GszeyJeGSub08yamNlsM1ttZqvMLCOkOUdHfuarzewxM7OAMvaPfJ8+MbN5Zta0yG1hev+UmDNk759Sv5eR28v3/nHO6aPYB1AV+AL4JVADWA60LHafDOAU4Hmga5Hj9YB1kcvDI9cPD2HO44HjItePAnKAw8KUscjtY4EpwBNh/JlHblsIXBC5XheoHbacwFnA4shzVAWWAh0Cynhe4fcI+BPwYuR62N4/peUM0/unxIxFbi/X+0cthZK1BdY659Y553YD04DLi97BObfBOfcJUFDssb8B5jjnvnXO/ReYA1wYtpzOuc+dc2si1zcDW4D9znRMdEYAM2sDNAJmxyFbTHKaWUugmnNuTuR+25xz28OWE3BALfwfl5pAdeA/AWVcUOR79C7QOHI9bO+fEnOG7P1T2veyQu8fFYWSHQ18VeTzTZFj8X5secXktcysLf4PxRcxylVUhTOaWRXgb8AdcchVXGW+l8cD35nZ383sIzMbY2ZVY57Qq3BO59xSYAH+v9oc4F/OudUxT1j+jH2AmRV8bGVUJmdUyN4/0YwVff9o57WSlXSetaxjdyvz2PKq9GuZ2ZHA/wP+4Jz72X/qMVCZjP2AfzrnvorTqe+iKpOzGtAeOBXYCLwIXAv8b0yS7avCOc2sOXAiP/0nOcfMznHOvRWrcIUvVcKxEjOaWS8gEzi3vI+NgcrkLDwemvdPCRkr9P5RUSjZJuCYIp83BjaX47Edij12YUxSlfxaFc2JmR0CvAkMdc69G+NshSqT8UygvZn1w5+nr2Fm25xzP+tsi4HK/sw/cs6tAzCzV4F2xKcoVCbnFcC7zrltAGY2E58z1kWhTBnN7H+AIcC5zrldRR7bodhjF8Y4X6HK5AzV+6eUjBV7/8S6YyQVPvDFch3QjJ86d04q5b6T+HlH83p8J9nhkev1QpizBjAPuC2s38tit11LfDuaK/O9rBq5f8PI588Cfw5hzt8CcyPPUT3y8780iIz4VtUXRDprixwP1ftnPzlD8/4pLWOx+7nh6nMAAAEXSURBVJT5/RO3LybZP4CLgc8j3+whkWP3AZdFrp+Or+I/AnnAyiKPvQ5YG/noHcacQC9gD/BxkY/WYcpY7DnK/Esd0M/8AuAT4FP8H+MaYcuJL14TgNXAKuDhADPOxXdyF/7uvV7ksWF6/5SYM2Tvn1K/l0Weo8zvHy1zISIiURp9JCIiUSoKIiISpaIgIiJRKgoiIhKloiAiIlGavCYSA2ZWHz9uHeAIYC+QG/m8rfPr1oiEnoakisSYmd0DbHPOPRR0FpHy0ukjERGJUlEQEZEoFQUREYlSURARkSgVBRERiVJREBGRKA1JFRGRKLUUREQkSkVBRESiVBRERCRKRUFERKJUFEREJEpFQUREolQUREQk6v8DrUH7PdaUS1QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#e over T\n",
    "ax.cla()\n",
    "data = pd.read_csv('./data/energy_T.csv')\n",
    "T_lqcd = data['T'] / 1000\n",
    "energy_lqcd = data['energy']\n",
    "ax.scatter(T_lqcd, energy_lqcd, color='cyan', label='LQCD')\n",
    "ax.plot(T, Ed/T**4, color='red', label='hrg')\n",
    "AddLabelAndLegend(ax, 'T', r'$\\epsilon/T^4$')\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEGCAYAAACJnEVTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAfS0lEQVR4nO3deXhU5fnG8e8joIiIC1D0B2KwakGxooTUSt0Qd+tSDSrohVSLLVbrgoqi1VqVKq11wyrWtaK0aqVWxIrIIi5oWFQE9wJSFiGogEEM4fn98U4gQAKTZOacyZz7c11zZfZzMzBPXt7zLubuiIhIMmwVdwAREYmOir6ISIKo6IuIJIiKvohIgqjoi4gkSOO4A2xOq1atvKCgIO4YIiINytSpU5e6e+vqHsvpol9QUEBJSUncMUREGhQzm1vTY+reERFJEBV9EZEEUdEXEUmQnO7Tr055eTnz58/n22+/jTtK5Jo2bUq7du1o0qRJ3FFEpIFqcEV//vz5bL/99hQUFGBmcceJjLtTWlrK/Pnz6dChQ9xxRKSBanDdO99++y0tW7ZMVMEHMDNatmyZyP/hiEjmNLiiDySu4FdK6p9bpCEYARQQimpB6nYuanDdOyIiuWYE0B8oS92em7oN0CeWRDVrkC39uM2ZM4fOnTvHHUNEcsRg1hf8SmWp+3ONin6WrFmzJu4IIhKRebW8P04q+nVUUVHBL37xC/bdd1+OPvpoVq1axeGHH84111zDYYcdxp133smnn37KQQcdRLdu3fjtb39L8+bN444tIlnQvpb3x6lh9+lfcgnMmJHZ9+zSBe64Y4tP+/jjj3nyySd54IEH6NWrF8888wwAX331FRMnTgTgxBNP5De/+Q1nnXUW9913X2ZzikjOuJkN+/QBmqXuzzVq6ddRhw4d6NKlCwBdu3Zlzpw5AJxxxhnrnvPGG29QXFwMQO/evSPPKCLR6AMMB3YHLPVzOLl3Ehcaeks/jRZ5tmyzzTbrrjdq1IhVq1YBsN1228UVSURi1IfcLPIbU0s/iw466KB13T4jR46MOY2IiIp+Vt1xxx3cfvvtFBUVsXDhQnbYYYe4I4lIwjXs7p2YFBQUMHPmzHW3Bw4cWO3z2rZty5tvvomZMXLkSAoLC6OKKCJSrciLvpntCPwV6Aw48HN3fyPqHFGYOnUqv/71r3F3dtxxRx566KG4I4lIwsXR0r8TeNHdTzezrQkjm/LSIYccwjvvvBN3DBGRdSIt+mbWAjgUOBfA3b8Dvosyg4hIkkV9IncPYAnwsJlNN7O/mtkGYxzNrL+ZlZhZyZIlSyKOJyKS36Iu+o2BA4G/uPsBwDfAoKpPcPfh7l7o7oWtW7eOOJ6ISH6LuujPB+a7+5TU7acJvwRERCQCkRZ9d18EfG5mP0jddSQwK8oMmVDTwmnDhw+nY8eOdOzYkcLCQiZMmLDusfLycgYNGsRee+1F586dKSoqYsyYMUAYArrffvux3377sc8++3DttdeyevXqKP4oIpIwcYzeuQgYkRq58xnQL4YMGff8889z//33M3nyZFq1asW0adM46aSTmDJlCm3btuW6665j4cKFzJw5k2222YbFixevW5gNYPz48bRq1YqVK1fSv39/+vfvz6OPPhrjn0hE8lHkM3LdfUaqz/6H7n6Ku3+ZzeNFtYXZrbfeytChQ2nVqhUABx54IP369WPYsGGUlZXxwAMPcPfdd69bs6dNmzb06tVrk/dp3rw59913H6NGjWLZsmVZSisiSZXXyzBUbmE2lzALrHILs2wU/vfff5+uXbtucF9hYSGzZs3ik08+oX379rRo0SKt92rRogUdOnTg448/zkJSEUmyvC76cW9h5u6xvFZEpCZ5XfSj3MJsn332YerUqRvcN23aNAoLC9lzzz2ZN28eK1asSOu9VqxYwZw5c9h7772zkFREkiyvi36UW5hdeeWVXHXVVZSWlgIwY8YMnn32WS644AKaNWvGeeedx8UXX8x334UJyAsXLuTxxx/f5H1WrlzJgAEDOOWUU9hpp52ykFREkiyvV9nM1hZmZWVltGvXbt3tyy67jMsuu4wFCxbQvXt31qxZw6JFi3jnnXeonGB20003ce2117LPPvvQtGlTtttuO2688cZ173HEEUfg7qxdu5ZTTz2V6667rp4pRUQ2Zbncd1xYWOglJSUb3Dd79mw6deqU9nuMIPThzyO08G8m+7vbrFmzhn79+rF27Voef/xxzCxj713bP7+IJI+ZTXX3atdyz+uWPsSzhVnjxo3529/+FvFRRUS2LK/79EVEZEMNsujncpdUNiX1zy0imdPgin7Tpk0pLS1NXAF0d0pLS2natGncUUSkAWtwffrt2rVj/vz5JHGt/aZNm24wakhEpLYaXNFv0qQJHTp0iDuGiEiD1OC6d0REpO5U9EVEEkRFX0QkQVT0RUQSREVfRCRBVPRFRBJERV9EJEFU9EVEEkRFX0QkQVT0RUQSREVfRCRBVPRFRBJERV9EJEEiX2XTzOYAK4AKYE1N+ziKiEjmxbW08hHuvjSmY4uIJJa6d0REEiSOou/AS2Y21cz6x3B8EZHEiqN7p7u7LzCz7wFjzewDd59U+WDqF0F/gPbt28cQT0Qkf0Xe0nf3BamfXwDPAkUbPT7c3QvdvbB169ZRxxORTJk1C55+Ou4UspFIi76ZbWdm21deB44GZkaZQUQisGYN9O0Lv/oVfP113Gmkiqhb+m2AyWb2DvAWMNrdX4w4g4hk29ChUFLCq/feS8EOO7AVUACMiDmWRNyn7+6fAftHeUwRidjMmXDDDcwtLubY4mLKUnfPJXWyDugTUzTRkE0RyaQ1a6BfP2jRglOHDVtX8CuVAYPjyCXrxDU5S0TyUapbh3/8gxk1DMSYF3Ek2ZBa+iKSGaluHYqLobiYmgZcayB2vFT0RaT+Krt1dtgBhg0D4Gag2UZPa5a6X+Kj7h0Rqb/bbgvdOk89BaluncqTtYMJXTrtCQVfJ3HjpaIvIvVTtVvn9NM3eKgPKvK5Rt07IlJ35eVw7rmw447runUkt6mlLyJ1N3QoTJ26QbeO5Da19EWkbiq7dXr12qRbR3KXir6I1F7Vbp177ok7jdSCundEpPZuu03dOg2UWvoiUjvvvQe/+526dRooFX0RSZ+6dRo8de+ISPpuugmmTVO3TgOmlr6IpGfixFD0+/ZVt04DpqIvIltWWgpnnw3f/z7cfXfcaaQe1L0jIpvnDuefD4sXwxtvwPbbx51I6kFFX0Q27/77YdQo+NOfoGvXuNNIPal7R0RqNnMmXHopHHMMXHJJ3GkkA1T0RaR6q1bBmWdCixbw6KOwlcpFPlD3johU7/LL4f334cUXoU2buNNIhuhXt4hs6tln4S9/gYEDQ9eO5A0VfRHZ0Oefw3nnhZO2N2tzw3yjoi8i61VUhPH45eXw5JOw9dZxJ5IMU5++iKx3yy0waVI4cbvXXnGnkSxQS19EgtdeC5ui9OkD55wTdxrJkliKvpk1MrPpZvZ8HMcXkY18+SX07g0FBXDvvWAWdyLJkri6d34DzAZaxHR8EankDhdcAAsWhNZ+C30t81nkLX0zawecAPw16mOLSDXuuisslXzTTVBUFHcaybI4unfuAK4E1lb3oJn1N7MSMytZsmRJtMlEkmbixDAJ6+ST4Yor4k4jEYi06JvZicAX7j61pue4+3B3L3T3wtbapEEke+bPD1sefv/78NhjWmYhIaLu0+8OnGRmxwNNgRZm9ri7nx1xDpFkW706bIRSVgYTJqgfP0Ei/dXu7le7ezt3LwDOBF5RwReJwUUXwZQpYTx+p05xp5EI6f9zIknzwAPhcvXV8LOfxZ1GIhbbjFx3nwBMiOv4Iok0ZQr8+tdhEbXf/z7uNBIDtfRFkmLxYjjtNGjbFp54Aho1ijuRxGCLLf3USVcAA84HHnD3F7KaSkQyq7wcioth2bKwz+3OO8edSGKSTkv/RmBfoBXQLPVTRBqSgQPh1Vfhr3+F/fePO43EKJ2ifyjQHPgWeN/dH8tuJBHJqMcfD7NuL7kkrK8jibbFou/uZe5+PbAUKMt+JBHJmOnT4Re/gMMOg9tuizuN5IC0Ru+Y2R/dfSAwLst5RCRTSkvDkMxWreAf/4AmTeJOJDkg3SGbPbKaQkQyq7wczjwzrJz56qvwve/FnUhyhHbOEsk37jBgALz8Mjz0kFbOlA2kW/T3N7P/Au8BM6v8/MDdy7MVTkTqYMiQMEpn8GDo1y/uNJJj0p2c9S5hsbR7gFLgaOBhYKmZzcxSNhGprSeeCMW+Tx/NuJVqpd294+4LgAXAS5X3mZkBe2Yhl4jU1qRJoWV/2GHw4IPa8lCqlW5L/9lUgd+ABx9nOJOI1NYHH8App8Aee8Czz8I228SdSHJUui39/wOmmtlHwIvAi+6+KHuxRCRtixfD8ceHIZkvvAA77RR3IslhaRV9d/8lgJl1BI4DHjGzHYDxhF8Cr7l7RdZSikj1ysrgpJNg0aKw9WGHDnEnkhxXq1U23f0Dd/+zux9LGLs/GSgGpmQjnIhsRkVFOGH79tvw5JPQrVvciaQBSKvom1lfM1tqZsvM7FEz297dV7n7C+5+kbsXZjuoiGzk8sth1Ci4886wsblIGtJt6V8HHAV0BOYBt2QtkYhs2Z13hssll4StD0XSlO6J3OXuPj11/TozU3eOSFxGjYJLL4VTT4U//jHuNNLApFv0dzWz/sBs4ANAKzeJxOGtt8LyyEVFYclk7X4ltZRu9871wA+Bm4APgc5m9oKZDTGzs7KWTiShRgAFhC9oQeo2M2fCccfBLrvAc89Bs2bxBZQGK90hm8Or3jazdoRfAvsBJwBPZj6aSDKNAPqzfvOKucBtH33Ez3r2ZNumTcNCalo1U+oo7SGbZlZkZpVjwloQTuq+5+5nZyWZSEINZsPdinafM4fnjzySb9auhXHjwqxbkTpKdxOV6wmTshqb2VjgR8AEYJCZHeDuN2cvokiyzKty/f/+9z9e6dGD7b75hh7jxzOjY8fYckl+SPdE7ulAF2AbYBHQzt2Xm9lQwsQsFX2RDGlP6NJp/cUXvNyzJ62WLuXIceP4ShuaSwak272zxt0r3L0M+NTdlwO4+ypgbdbSiSTQzUDbZcsYe9RR7D53LieMHs2sbt3UspKMSLfof2dmlUMFulbemVp/R0VfJIP6LF/OjGOPpeMHH3DKv/7F54ccwnCgT9zBJC+k271zqLuvBnD3qkW+CdA33YOZWVNgEqGbqDHwtLtfn+7rRfLeN9/ACSfQavp0+Oc/eemoo+JOJHkm3SGbq2u4fymwtBbHWw30cPeVZtYEmGxmY9z9zVq8h0h++vbbsCb+66+HBdR++tO4E0keinRjdHd3YGXqZpPUxaPMIJKTysuhV68wBv+RR8J1kSyo1dLKmWBmjcxsBvAFMNbdp2z0eH8zKzGzkiVLlkQdTyR6a9aEJZL//W8YNgz6pt1jKlJrkRf91CigLkA7oMjMOm/0+HB3L3T3wtatW0cdTyRa5eVwzjnw1FMwdCgMGBB3IslzkRf9Su7+FWGC17FxZRCJ1apV8LOfwciR8Ic/wMCBcSeSBIi06JtZazPbMXV9W6AnYdVOkWRZsQJOOAFGj4Z774Wrroo7kSREpCdygV2BR82sEeEXzj/c/fmIM4jEa9mysJF5SQk89hicreWrJDpRj955FzggymOK5JTFi+Hoo+GDD+CZZ7TNoUQu6pa+SHLNmwc9e8L//he6dXr2jDuRJJCKvkgUPv4YjjwSli+HsWPh4IPjTiQJpaIvkm3vvhu6dCoqYPx4OEA9nBKf2IZsiuSbarc4nDIFDjsMGjeGV19VwZfYqeiLZEDlFodzCeuKzAUeHz+e8iOPhJYtYfJk0AYokgNU9EUyYOMtDk97+mmePe44PisoCC38goJ4golsREVfJAPWbXHozqAhQ3i6uJipXbvSfeJE2HXXOKOJbEBFXyQD2gNNvvuOh/v1Y8g11zCid2+OHDeO5i1bxh1NZAMavSOSAUNLS2lz2mkcOnEi199wAzf+9rc0M9MWh5JzVPRF6uujjyg+8UQq5s7l4hEjuKd3b3Yn7HWrLQ4l16joi9THhAlhpcxGjWj0yivc1b07d8WdSWQz1KcvUlcPPxwmXe2ySxiP37173IlEtkhFX6S21q6Fa66Bn/88TLx6/XXYY4+4U4mkRd07IrVRVha2M3z6abjgArj7bmjSJO5UImlT0RdJ17x5cNppMHUq/OlPcOmlYBZ3KpFaUdEXSceYMWGzk/JyGDUKTjop7kQidaI+fZHNqaiAa68NO13tthtMm6aCLw2aWvoiNVm0CHr3Dsshn3de6L/fdtu4U4nUi4q+SHUmToQzz4Svv4ZHHgknb0XygLp3RKpauxZuvRV69IAWLcL4exV8ySMq+iIpTy1bxssnnwyDBvHv4mL+XlIC++0XdyyRjFL3jggw5u23KSouZtcFC7jwnnu4d8AAmpmxBq2fI/lFLX1JtrVr4e676fGTn+DATyZP5t4LLwQzygibo4jkExV9Sa7PP4djjoGLL+alo4/mwGnTeLuoaIOnzKvhpSINlYq+JI97GJHTuTO88Qbcfz8XPfccX+688yZPbR99OpGsUtGXZFm0CE45Bfr1gy5d4N13oX9/bjaj2UZPbQbaBEXyTqRF38x2M7PxZjbbzN43s99EeXxJuKeeCq37//wHbr89TLpKrY7ZBxgO7A5Y6udwdBJX8k/Uo3fWAJe7+zQz2x6YamZj3X1WxDkkSZYtgwsvhJEjoVs3ePRR6NRpk6f1QUVe8l+kLX13X+ju01LXVwCzgbZRZpCEGT0a9t03LIX8+9+Hte+rKfgiSRFbn76ZFQAHAFM2ur+/mZWYWcmSJUviiCb5YPlyOP98OPFEaN0a3n47LJzWWFNTJNliKfpm1hx4BrjE3ZdXfczdh7t7obsXtm7dOo540pC5w8iRlHXsSMXDDzPk6qvZ++23GdGlS9zJRHJC5M0eM2tCKPgj3P2fUR9f8tiHH4a++3Hj+PDAA7lg1Kh14+77p56iPntJuqhH7xjwIDDb3W+P8tiSx8rKwp61++0HJSVcO2wYhW+9tcFEK82uFQmi7t7pDpwD9DCzGanL8RFnkHzhHnax6tQJhgyBs86CDz/klgEDWNuo0SZP1+xakYi7d9x9MmEYtEj9fPYZXHxxGJ3TuTNMmgSHHAKEWbRzq3mJZteKaEauNDTffgs33hiGYU6cGDYonzZtXcGHMItWs2tFqqfxa9IwuMO//gVXXAGffAJnnBEKfttNp3lUnqwdTOjSaU8o+DqJK6KiLw3Bq6/CVVeFxdE6doSxY6Fnz82+RLNrRaqn7h3JXTNnwk9/CoceStncuVz1wAM0ee89Cnr2ZETc2UQaKLX0JffMmwfXXx/WyGnRgulDhnDUxRdT2iz01M9F4+5F6kotfckdpaUwcCDsvTc8+SRcfjl89hmnDhq0ruBX0rh7kbpRS1/iV1YGd94Jt94a1szp2xd+9ztoHwZZ1jS+XuPuRWpPRV/i8803MHw4DB0KCxeG/vtbbgnj7qvQuHuRzFH3jkTvyy/DMse77w6XXRa6cyZNguee26Tgg8bdi2SSWvoSnUWL4M9/hnvvhZUrw7LHV18NBx+82Zdp3L1I5qjoS/bNmRO6cB58EMrLoVcvGDQI9t8/7bfQuHuRzFDRl+yZPRv+8AcYMQK22grOPReuvBL23DPuZCKJpaIvmeUOEybAXXeFZRO23RYuuigMv2zXLu50IomnE7mSGStWhL76zp2hR4+wdMLgwaFr589/hnbtGAEUEP7RFYBm1YrEQC19qZ8PP4Rhw+CRR0LhP/BAePjhsCDattuue9oIwizastRtzaoViYeKvtReRUVYx/6ee8LiZ02ahCJ/4YXwox+BbbplwmDWF/xKlbNqVfRFoqOiL+lbuhQeeih048ydG5Y1vukmOP98aNNmsy/VrFqR3KCiL5v33XcwZkzovhk9Ogy5PPzwsJb9ySdD4/T+CWlWrUhuUNGXTbnD9Olhlcsnnggt/DZtwiicfv2qnTW7JTezYZ8+aFatSBxU9GW9RYvCmPpHH4X33oOttw6t+b594Zhj0m7VV0ezakVyg4p+0q1aBc8/Hwr9iy+Gk7Q/+lHotz/jDNh554wdSrNqReKnop9EZWWhwD/1VCj4K1eGk7JXXBFa9R07xp1QRLJERT8pvvkGXnghFPrRo0Phb9UKeveG4mI44gho1AgIY+rVDSOSn1T089mKFaHAP/VUGIGzalU4Idu3L5x+Ohx66Cb99JpEJZLfVPTzzcKFoUX/73+HLpzVq2HXXeG880Kh/8lP1rXoq6NJVCL5TUW/oauogLffDoV+9GiYNi3cv9tu8MtfhkJ/8MFhlcs0aBKVSH6LtOib2UPAicAX7l77wd4SfPkl/Oc/odCPGRPG0W+1VSjuQ4bACSeEsfTVLIewJZpEJZLfom7pPwLcAzwW8XEbtoqK0IIfNy4U+tdfD/e1bAnHHQfHHx/G0WdgeKUmUYnkt0iLvrtPMrOCKI/ZILnDrFnwyiuh0E+YAF9/HR7r0iXsOnXCCVBUtNn++brQJCqR/JZzffpm1p/UgJH27RPUqTBnTijw48aFYr94cbi/Q4cwpLJHj3DZwsJmmaBJVCL5K+eKvrsPB4YDFBYWesxxssMdPvkEJk8Ol/Hj4b//DY+1aQNHHrm+yHfoEG9WEckrOVf081J5eVjArLLIv/YafPFFeGznneGQQ+DSS0Ox79RpkxOwmiwlIpmiop8NX30FU6asL/JTpoSJUQB77AHHHhvGy3fvHpY82MxwSk2WEpFMinrI5pPA4UArM5sPXO/uD0aZIeNWr4Z33oG33lp/+fDD8FijRuHEa//+64v8rrvW6u01WUpEMinq0TtnRXm8jFu7Fj76aMMCP2NG6L4B2GWXsELlOeeEnwcdBM2b1+uQmiwlIpmk7p10PfZY2ERk+fJwu3lz6NYNLrssDJ0sKgorVdZhQtTmaLKUiGSSin669toL+vRZX+B/8IOMj5GvjiZLiUgmqehvxgajZn78Y27+8Y8j70fXZCkRySQV/Rrk0qgZTZYSkUxJb+nFBNrcqBkRkYZKRb8GGjUjIvlIRb8GNY2O0agZEWnIVPRrcDNhlExVGjUjIg2din4N+hBWfdsdsNTP4eiEqog0bBq9sxkaNSMi+UYtfRGRBFHRFxFJEBV9EZEEUdEXEUkQFX0RkQRR0RcRSRBzz929x81sCdUvJx+XVsDSuENsgTJmhjJmhjJmRm0z7u7urat7IKeLfq4xsxJ3L4w7x+YoY2YoY2YoY2ZkMqO6d0REEkRFX0QkQVT0a2d43AHSoIyZoYyZoYyZkbGM6tMXEUkQtfRFRBJERV9EJEFU9FPM7Fgz+9DMPjGzQdU8fqiZTTOzNWZ2+kaP9TWzj1OXvrmW0cy6mNkbZva+mb1rZmfkWsYqj7cws/+Z2T25mNHM2pvZS2Y228xmmVlBDma8LfV3PdvM7jIziynjZanP6F0zG2dmu1d5LFe+M9VmzLHvTI2fY+rx2n1n3D3xF6AR8CmwB7A18A6wz0bPKQB+CDwGnF7l/p2Bz1I/d0pd3ynHMu4N7JW6/n/AQmDHXMpY5fE7gSeAe3Lt7zr12ATgqNT15kCzXMoIHAy8lnqPRsAbwOExZTyi8vMBfgX8PXU9l74zNWXMpe9MtRmrPF6r74xa+kER8Im7f+bu3wEjgZOrPsHd57j7u8DajV57DDDW3Ze5+5fAWODYXMro7h+5+8ep6wuAL4BqZ+vFlRHAzLoCbYCXspCt3hnNbB+gsbuPTT1vpbuX5VJGwIGmhAKyDdAEWBxTxvFVPp83gXap67n0nak2Y459Z2r6HOv0nVHRD9oCn1e5PT91X7ZfWxsZOY6ZFREKwqcZylVVnTOa2VbAn4ArspCrqvp8jnsDX5nZP81supkNNbNGGU9Yj4zu/gYwntAyXQj8x91nZzxh7TOeB4yp42vrqj4Z18mx78y6jHX9zmi7xKC6Ps90x7LW57W1Ue/jmNmuwN+Avu6+SUs7A+qTcQDwgrt/nqUu6Er1ydgYOAQ4AJgH/B04F3gwI8nWq3NGM9sT6MT61uBYMzvU3SdlKlzloaq5r9qMZnY2UAgcVtvX1lN9MlbenzPfmWoy1uk7o6IfzAd2q3K7HbCgFq89fKPXTshIqk2PU9eMmFkLYDRwrbu/meFsleqT8cfAIWY2gNBXvrWZrXT3TU5sxZhxPjDd3T8DMLNRwEFkvujXJ+OpwJvuvhLAzMYQMma66KeV0cx6AoOBw9x9dZXXHr7RaydkOF99M+bUd6aGjHX7zmT6xERDvBB++X0GdGD9yZR9a3juI2x6Ive/hBNSO6Wu75xjGbcGxgGX5OrnuNFj55K9E7n1+RwbpZ7fOnX7YeDCHMt4BvBy6j2apP7efxpHRsL/iD4ldUK0yv05853ZTMac+c7UlHGj56T9ncnaH6ahXYDjgY9SH+7g1H03Aielrncj/Fb+BigF3q/y2p8Dn6Qu/XItI3A2UA7MqHLpkksZN3qPtP8Bx/B3fRTwLvAeoeBunUsZCb+Y7gdmA7OA22P8HF8mnESu/Df3XA5+Z6rNmGPfmRo/x7p8Z7QMg4hIgmj0johIgqjoi4gkiIq+iEiCqOiLiCSIir6ISIJocpZILZhZS8L4bYBdgApgSep2kYf1U0RyloZsitSRmd0ArHT3P8adRSRd6t4REUkQFX0RkQRR0RcRSRAVfRGRBFHRFxFJEBV9EZEE0ZBNEZEEUUtfRCRBVPRFRBJERV9EJEFU9EVEEkRFX0QkQVT0RUQSREVfRCRB/h+a/qsCfnPjIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P over T\n",
    "ax.cla()\n",
    "data = pd.read_csv('./data/pressure_T.csv')\n",
    "T_lqcd = data['T'] / 1000\n",
    "pressure_lqcd = data['pressure']\n",
    "ax.scatter(T_lqcd, pressure_lqcd, color='cyan', label='LQCD')\n",
    "ax.plot(T, 3*Pr/T**4, color='red', label='hrg')\n",
    "AddLabelAndLegend(ax, 'T', r'$3P/T^4$')\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3zV9fXH8dcxgIjKRhAiU6zgQonUShWxakUUcDEEB1KxVou71Z8odbWOCuICYx0gCNUiGtwDwUkliIOhMowSmYIiyAyc3x/fm/SS3NxcktyVvJ+Px33k3u+655txTz7b3B0REZHi9kh2ACIikpqUIEREJCIlCBERiUgJQkREIlKCEBGRiGokO4DK1LhxY2/dunWywxARSRtz5sz5wd2bRNpXpRJE69atyc3NTXYYIiJpw8y+LW2fqphERCQiJQgREYlICUJERCKKWxuEmT0BnA6sdvdDI+y/HhgYFkcHoIm7rzOzPGADsAMocPes8saxfft28vPz2bJlS3kvkZZq165NZmYmNWvWTHYoIpKm4tlI/RTwEDA+0k53vxe4F8DMzgCudvd1YYd0d/cfKhpEfn4+++67L61bt8bMKnq5tODurF27lvz8fNq0aZPscEQkTcWtisnd3wXWlXlgYAAwKR5xbNmyhUaNGlWb5ABgZjRq1KjalZpEpHIlvQ3CzOoApwJTwjY78IaZzTGzoWWcP9TMcs0sd82aNaUdU2nxpovqeM8iUrmSniCAM4APilUvdXX3o4AewOVmdnxpJ7t7trtnuXtWkyYRx3qIiFRdr70GDzwA27ZV+qVTIUH0p1j1krsvD31dDUwFuiQhrkqRl5fHoYeWaKMXEak4dxgxAh58EDIyKv3ySU0QZlYP6Aa8GLZtbzPbt/A5cAowLzkRJkZBQUGyQxCRdPTBB/Dxx3D11XFJEPHs5joJOAFobGb5wAigJoC7jw0ddibwhrv/EnZqU2BqqA69BvCMu78WrzgTYceOHVxyySV8+OGHtGjRghdffJEePXpw7LHH8sEHH9CrVy/69OnDwIED2bFjBz169GDkyJFs3Lgx2aGLSCr75z+hUSO46KK4XD5uCcLdB8RwzFME3WHDty0FjohLUFddBZ9+WrnX7NQJ7r8/6iGLFi1i0qRJPPbYY/Tt25cpU4L2+J9++omZM2cCcPrpp3PllVcyYMAAxo4dG+1yIiKwaBHk5MBNN0GdOnF5i1Rog6jy2rRpQ6dOnQDo3LkzeXl5APTr16/omI8++ohzzz0XgPPOOy/hMYpImhk1CmrWhMsvj9tbVKnZXMtUxn/68bLnnnsWPc/IyGDz5s0A7L333kmJR0TS3A8/wFNPwaBB0KxZ3N5GJYgUccwxxxRVPU2ePDnJ0YhIShszBjZvhmuuievbKEGkiPvvv5+RI0fSpUsXVqxYQb169ZIdkoikoi1b4KGHoEcPOOSQuL5V9apiSoLWrVszb97/euled911EY9r0aIFs2bNwsyYPHkyWVnlnp9QRKqyiRNh9Wq49tq4v5USRIqYM2cOV1xxBe5O/fr1eeKJJ5Idkoikmp074b774Igj4MQT4/52ShAp4rjjjuOzzz5Ldhgikspeew0WLoSnn4YEzLemNggRkXRx333QogWEdZGPJyUIEZF0MHcuTJ8Ow4YF4x8SQAlCRCQd3Hcf7LMPDI26AkKlUoIQEUl1+fnw73/DH/4A9esn7G2VIBJgn332ibg9Ozubgw8+mIMPPpisrCxmzJhRtG/79u3ccMMNtG/fnkMPPZQuXbrw6quvAkHX2cMOO4zDDjuMjh07Mnz4cLZu3ZqIWxGRZHjggaAH05VXJvRtlSCS5KWXXuLRRx/l/fff58svvyQ7O5tBgwbx/fffA3DzzTezYsUK5s2bx7x585g2bRobNmwoOv+dd97hiy++4OOPP2bp0qUMTWCxU0QS6Oef4dFH4ZxzoHXrhL61EkQxE4HWBN+Y1qHX8XD33Xdz77330rhxYwCOOuooBg8ezMMPP8ymTZt47LHHePDBB4vmcWratCl9+/YtcZ199tmHsWPH8sILL7BuXaxLgItI2nj88SBJJGBgXHFKEGEmAkOBbwkWxf429DoeSWL+/Pl07tx5l21ZWVksWLCAxYsX07JlS+rWrRvTterWrUubNm1YtGhRHCIVkaQpKIDRo+G446BL4hfWVIIIcxOwqdi2TaHtieDuSTlXRFLUlCnw7bdJKT2AEsQuvtvN7RXRsWNH5syZs8u2Tz75hKysLA488EC+++67XdocotmwYQN5eXkcdNBBcYhURJLCPVgxrn17OOOMpISgBBGm5W5ur4i//OUv/PWvf2Xt2rUAfPrpp0ydOpVLL72UOnXqMGTIEIYNG8a2bdsAWLFiBRMmTChxnY0bN/KnP/2JPn360KBBgzhEKiJJ8d57kJsbrDe9R3I+qjUXU5g7CdocwquZ6oS2V8SmTZvIzMwsen3NNddwzTXXsHz5crp27UpBQQErV67ks88+o0mTJgDccccdDB8+nI4dO1K7dm323ntvbrvttqJrdO/eHXdn586dnHnmmdx8880VjFJEUsp99wXrTV94YdJCsKpUd52VleW5ubm7bFu4cCEdOnSI+RoTCdocviMoOdwJDKzEGCMpKChg8ODB7Ny5kwkTJmCVNAnX7t67iKSIr76CDh1g+HAI+8cwHsxsjrtHXF8gbiUIM3sCOB1Y7e6HRth/AvAi8E1o0/Puflto36nAaCAD+Je73xWvOIsbSPwTQnE1atTg6aefTvC7ikjKGjUKatWK63rTsYhnxdZTwKllHPOeu3cKPQqTQwbwMNAD6AgMMLOOcYxTRCR1rFkD48YF6003bZrUUOKWINz9XaA8I7e6AIvdfam7bwMmA70rGEtFTk9L1fGeRaqEMWOCZUXjvN50LJLdi+k3ZvaZmb1qZoWLq7YAloUdkx/aFpGZDTWzXDPLXbNmTYn9tWvXZu3atdXqA9PdWbt2LbVr1052KCKyOzZvDtabPu006Jj8ipNk9mL6BGjl7hvN7DTgBaA9EKmFttRPd3fPBrIhaKQuvj8zM5P8/HwiJY+qrHbt2rv0nBKRNDBhQlDFlKSBccUlLUG4+89hz18xs0fMrDFBieGAsEMzgeXlfZ+aNWvSpk2b8gcqIpIIO3fCyJHQqRN0757saIAkJggzawascnc3sy4E1V1rgZ+A9mbWBvge6A+cl6w4RUQS4tVX4csvg1JEAtabjkU8u7lOAk4AGptZPjACqAng7mOBc4DLzKwA2Az096ChoMDMrgBeJ+jm+oS7z49XnCIiKeGf/4TMTIgwa3OyxC1BuPuAMvY/BDxUyr5XgFfiEZeISMrJzYUZM+CeexK23nQskt2LSUREbr4ZGjRI6HrTsdBcTCIiyfTuu/Daa3D33VCvXrKj2YVKECIiyeION94I++8PV1yR7GhKUAlCRCRZXn4ZPvwwGD1dp06yoylBJQgRkWTYuRP+7/+gXTsYMiTZ0USkEoSISDJMmgRffAHPPJNSPZfCqQQhIpJo27bBLbfAEUdAv37JjqZUShAiIpVkItCa4IO1deh1RI8/DkuXwp13Jm050VioiklEpBJMZNcli78NvYZii5Bt2gS33w5duwaztqaw1E1dIiJp5Ep2Xc+e0Oubih/44IOwYgXcdVfKzLlUGiUIEZEKmkgw02gk34W/+PHHIDGcdhr89rfxD6yClCBERCqoRCkhTMvwF/feCz/9FLQ9pAElCBGRCvo2yr6iVLByJYweDf37B2s+pAElCBGRCphI5GUwARoR1kB9xx1B99bbb09IXJVBCUJEpAKuJPKayAaMLnzxzTeQnR2MmD7wwITFVlFKECIi5RStcdoJKz2MGAEZGcG03mlECUJEpBwmAhdG2d+q8MkXXwTLiP75z9CiRfwDq0RKECIiu6lwUNyOKMcUNU4PHw5168INN8Q9rsqmBCEispsiDYoLV9Q4/dFHkJMD118PDRsmJLbKpAQhIrIborU7ANQh1DjtHkznvd9+cOWVCYmtssUtQZjZE2a22szmlbJ/oJl9Hnp8aGZHhO3LM7MvzOxTM8uNV4wiIrsr2qC4DCCbUOnhzTdhxoygimmffRIRWqWLZwniKeDUKPu/Abq5++HA7QTf13Dd3b2Tu2fFKT4Rkd0WbVDcOELJoXAxoFatYOjQKGektrjN5uru75pZ6yj7Pwx7OQvIjFcsIiKVoXBQXKRxD7sMipsyBebMgXHjYM89ExVepUuVNoghwKthrx14w8zmmFnU9GtmQ80s18xy16xZE9cgRaR6i2lQXEFBMN7hkENg4MAIR6ePpK8HYWbdCRJE+NSGXd19uZntB7xpZl+6+7uRznf3bELVU1lZWZF+diIiFRbzoLhx4+Crr2Dq1GBwXBpLagnCzA4H/gX0dvei7727Lw99XQ1MBbokJ0IRkUC0xumiQXFbtsDf/ga//jX07h3/oOIsaQnCzFoCzwPnu/vXYdv3NrN9C58DpwARe0KJiCTCRGKcsXXMGMjPh3/8I+UXA4pF3KqYzGwScALQ2MzygRFATQB3HwvcQtCu84gF38iCUI+lpsDU0LYawDPu/lq84hQRiaZw1HRpihqnf/wR/v53OPlk6N49IbHFWzx7MQ0oY/8fgD9E2L4UOKLkGSIiiRdt1HTRoDiAG28MksQ99yQkrkQos4rJzE43s1Tp7SQikjBljZouGhT34Yfw6KPBiOk0WQwoFrF88PcHFpnZPWbWId4BiYikimgTZLQilBy2b4dLL4UDDoBbb01MYAlSZhWTuw8ys7rAAOBJM3PgSWCSu2+Id4AiIslQVumhqGF61CiYNw9efDFtp9QoTUxVR+7+MzAFmAzsD5wJfGJmf45jbCIiSROtW2tRw3ReXtCttU8f6NUrEWElVCxtEL3MbCownaAXUhd370HQkHxdnOMTEUmKaN1ai2ZrvfzyYDDcAw8kKKrEiqUX0znAqOIjmd19k5ldHJ+wRESSJ6Y5l6ZMgVdegZEjg/aHKiiWKqYVxZODmd0N4O5vxyUqEZEkKVxKNOqcS+vXw7BhcOSRwVKiVVQsCeLkCNt6VHYgIiLJVtZSokVzLg0fDitXBl1bayR9Sru4KfXOzOwy4E9AOzP7PGzXvsAH8Q5MRCTRylpKtBXA7Nnw8MNB+8PRRycmsCSJlvqeIZiC+x9A+GrbG9x9XVyjEhFJsFiWEv17QUEw5mH//eGOOxIUWfJESxDu7nlmdnnxHWbWUElCRKqSWJYSPe/BB2HuXHjuOahXL0GRJU9ZJYjTgTkEVW/hUxM60DaOcYmIJExZs7WOAwYuWxYsBHTaaXD22QmKLLlKTRDufnroa5vEhSMiklgxz9Y6bFiw1vTDD1eJqbxjEctAua6hdRkws0FmNjK0loOISNqLabbWF1+EF14IRk23bp2gyJIvlm6uY4BNZnYE8BeCktjTcY1KRCQBYpqtdePGYKzDYYfB1VcnKLLUEEuCKHB3B3oDo919NEFXVxGRtFU4IK40RbO1jhgBy5bB2LFQs2ZCYksVsYzw2GBmNwKDgOPNLIPQynAiIumorAFxEJqtde5cGD066Np67LEJiS2VxFKC6AdsBYa4+0qgBXBvXKMSEYmjsgbENQIG7tgRJIZGjYI1pquhWNaDWAmMDHv9HTA+nkGJiMRLLAPiRkNQpTR7NkycCA0aJCS2VBNLL6azzGyRma03s5/NbIOZ/ZyI4EREKlssA+IGLl8O//d/cNJJMGBAgiJLPbFUMd0D9HL3eu5e1933dfe6sVzczJ4ws9VmNq+U/WZmD5jZYjP73MyOCtt3YSgxLTKzaG1JIiIxiWlAHMBVV8HWrTBmTLUZ8xBJLAlilbsvLOf1nwJOjbK/B9A+9BhK0KUWM2sIjAB+DXQBRphZ9SzjiUiliHlA3KuvBlNpDB8OBx6YkNhSVSy9mHLN7N/ACwSN1QC4+/Nlneju75pZ6yiH9AbGh7rRzjKz+ma2P3AC8GbhfE9m9iZBopkUQ7wiIiXENCBu3bqgYfrgg+H66xMWW6qKJUHUJfi+nhK2zYEyE0QMWgDLwl7nh7aVtr0EMxtK6B+Dli01wFtESoppQJw7DB4crPPw4Yew554Jii51xdKLaXAc3z9S5V7xiQHDt5fc6J5N8PMlKysr4jEiUn3FPCDuoYcgJydYQjQrKyGxpbpYejEdZGZvFzY0m9nhZja8kt4/HwhfzDUTWB5lu4hIzGIeEPfJJ3DddXD66UEDtQCxNVI/BtwIbAdw98+B/pX0/jnABaHeTMcA6919BfA6cIqZNQg1Tp8S2iYiErObiGFA3IYN0K8fNGkCTz5ZrXstFRdLG0Qdd//Ydv2mFcRycTObRNDg3NjM8gl6JtUEcPexwCvAacBigp/j4NC+dWZ2OzA7dKnbtECRiOyuaF1a6wCj3eGyy2DpUnjnHWjcOFGhpYVYEsQPZtaOUBuAmZ0DrIjl4u4edYRJqPdSiRXrQvueAJ6I5X1ERIqbSNCYGalhsmhA3LhxwUjpW2+F449PaHzpIJYEcTnB9/JgM/se+IZQm46ISCoqbJiOlByM0IC4hQvh8suhe3e4Kdr46uorll5MS4GTQosG7eHuG+IflohI+ZTVMO3AwM2boX9/qFMHJkyAjIzEBZhGoiYIM/sVwff64NCmhWaW7e5fxz0yEZFyKKthuhXAtdfC55/DK69A8+aJCSwNldqLycx+A8wANhBUMT0G/ALMCPU4EhFJOWU1TI+fMiWYY+m666BHj0SFlZailSBuAQa4+4ywbS+Y2XSC3kj6zopIyphIMJ1GaTKAiXl5HD9kCHTpAnfemaDI0le0cRDtiiUHANx9JtA2bhGJiOymwnaH0qbTMGD89u30GTAA3GHSJKhVK3EBpqloJYhojdG/VHYgIiLlVVa7gwPn3XwzzJoFzz4LbfU/biyiJYgDzOyBCNuNUibOExFJtLLWeAA4//XX4e67g5lazz03EWFVCdESRLS5bnMrOxARkd1V1hoPAG1XrCD7/PPh0ENh1KhEhFVllJog3H1cIgMREdkdhYPhok3E12THDj4YNIjaGzfCjBmw116JCa6KiGWyPhGRlBLLLK0TgNV33UWz6dODqbw7dkxMcFWIEoSIpJ1oq8NBaI2H99+HW26BAQOChYBktylBiEhaKWt1uDrAvevWBYmhTRsYO1ZTeJdTmXMxmVkT4BKgdfjx7n5x/MISESmprNXhMoB/bdvGuX37wqpV8NFHULdugqKremKZzfVF4D3gLaJX+YmIxE0s7Q7j3BlwySXw9tswbhx07pyg6KqmWBcM+mvcIxERiSKm1eFGjIDx4+G22+CCCxIUWdUVSxvES2Z2WtwjEREpRVmD4eoALz3+ONx+OwwZAsOHJyiyqi2WBHElQZLYYmYbQo+f4x2YiAiUPRguA5j22mscc+ml8PvfBzO1qlG6UsSyYNC+iQhERKS4sgbD1QGemzuXE889Fw47DJ57DmrWTFyAVVwsbRCYWS+gcMHWGe7+UvxCEhGJcTDcd99xWs+e0KABvPwy7Kv/ZytTmVVMZnYXQTXTgtDjytC2MpnZqWb2lZktNrMbIuwfZWafhh5fm9lPYft2hO3Lif2WRCTdFZYcojVKH/bTT5zZowds2gSvvqqV4eIglhLEaUAnd98JYGbjgLlAiQ/8cGaWATwMnAzkA7PNLMfdFxQe4+5Xhx3/Z+DIsEtsdvdOsd6IiFQNsZQc6m/dyhtnngmLFsHrr8MhhyQouuol1pHU9cOe14vxnC7AYndf6u7bgMlA7yjHDwAmxXhtEamiyppGI8OdT4YModmMGfDkk9C9e4Iiq35iSRD/AOaa2VOh0sMc4O8xnNcCWBb2Op9S1pEws1ZAG2B62ObaZpZrZrPMrE9pb2JmQ0PH5a5ZsyaGsEQkVcUyjcanw4fTZuLEYMnQgQMTFFn1FLWKycwMeB84BjiaYLGgv7r7yhiuHamfmZdybH/gP+4eXqps6e7LzawtMN3MvnD3JSUu6J4NZANkZWWVdn0RSXGxTKMxPTubQ//+d7jkErjxxgRFVn1FTRDu7mb2grt3Bna3oTgfOCDsdSawvJRj+wOXF3vv5aGvS81sBkH7RIkEISLpL5Z2hzdfeYVf/+lPcNpp8MgjGuuQALFUMc0ys6PLce3ZQHsza2NmtQiSQIkkY2a/AhoAH4Vta2Bme4aeNwa6EvSgEpEqqKxpNE6cM4fuffvCEUfAv/8NNWLqoS8VFEuC6A58ZGZLzOxzM/vCzD4v6yR3LwCuAF4HFgLPuvt8M7stNK6i0ABgsruHVw91AHLN7DPgHeCu8N5PIlI1TCSYJjraNBoH5+XxUs+e0LgxvPQS7LNPYoITbNfP5QgHBA3IJbh7WeuEJ1xWVpbn5mq5bJF0UFitFHUCvh9/ZMmxx1Jv5Ur44AOtChcHZjbH3bMi7Su1nGZmDUNPN8QlKhGp1srqztpgyxbm9+lDvaVL4Y03lBySIFpF3hyCXkcGtAR+DD2vD3xH0C1VRGS3ldWd9eBffmF67940ffddmDgRunVLVGgSptQE4e5tAMxsLJDj7q+EXvcATkpMeCJS1ZTVnfXQ9ev5omfPYDW4cePgvPMSFZoUE0sj9dGFyQHA3V8FlM5FZLeV1Z214dq1zPjd7+C//w16K2nRn6SKpa/YD2Y2HJhAUOU0iOilQxGREsqaurvpypVMP/lkGi1aBC+8AD17JjA6iSSWEsQAoAkwFXgB2C+0TUQkJmWVHDKXLeO944+n/TffwCuvKDmkiFgWDFpH0OFARGS3lVVyaLtkCW//7nc0/+knar7xBhx7bAKjk2jKTBBm1gT4C3AIULtwu7ufGMe4RKQKKKvk0GHBAt4+6SQabNtGrenT4aijEhidlCWWKqaJwJcE3VpvBfIIptEQEYmocIT0IEof63DEp58ys1s36rtTe+ZMJYcUFEuCaOTujwPb3X2mu19MMLuriEgJhaWGaFMt/HrWLN7p3p06deqw13vvacGfFBVLL6btoa8rzKwnwYysmfELSUTSVVntDQDdZsxg2hln4M2asfdbb0GriLP5SAqIJUHcYWb1gGuBB4G6wNXRTxGR6iaWKbt//9prTD3zTLa2bUv9t96C/fdPUHRSHrH0Ynop9HQ9wcyuIiIllDVl95nPP8/k/v3ZeOihNHzjjWB2Vklp0Sbre5DSV4DD3YfFJSIRSTsTid7mcPGECTx20UXs0aULDV95BerXj3K0pIpoJYjwebNvBUbEORYRSUOFVUuluTQ7mzF//CN2wgmQk6P1HNJItMn6xhU+N7Orwl+LiED0Ruka27dz//XXc/no0cEyof/5D+y1V4IjlIqIdd2+6KsKiUi1E61Rusnq1Tzbty8nzJwJV10F99wDNWsmOEKpKC3sKiK7ZSJBg3RpbQ5Zs2fz/Fln0eSHH2DCBBg4MIHRSWUqdaCcmW0ws5/N7Gfg8MLnhdsTGKOIpICJQGOC0dGlJYeLnnyS9447jp0ZGUz/8EMlhzQXrQ1i30QGIiKpq6z1o2tu28aoq6/m8kce4a2TTuKnSZM4R91Y014sU22Um5mdamZfmdliM7shwv6LzGyNmX0aevwhbN+FZrYo9Ii2AJWIxFFhQ3RpyaHZihVMP/FELn/kEUZdfz2rX31VyaGKiFsbhJllAA8DJwP5wGwzy3H3BcUO/be7X1Hs3IYE3WqzCBrI54TO/TFe8YpISWWNjj7mo4+YcvbZ1Fu/nismT+Y3/fqhBUKrjniWILoAi919qbtvAyYDvWM89/fAm+6+LpQU3gROjVOcIlJMLLOxXpKdzcxu3diy117MmDWLh/r1Qy0OVUs8E0QLYFnY6/zQtuLONrPPzew/ZnbAbp6LmQ01s1wzy12zZk1lxC1SrZU1G2utrVvJvuQSsi+9lJknncSc3Fx6HnZYAiOURIlngrAI24qPp5gGtHb3w4G3gMLBeLGcG2x0z3b3LHfPatKkSbmDFanuYik1tMjPZ2a3blzyr38x76abOHnaNM5t0CBhMUpixTNB5AMHhL3OJJgqvIi7r3X3raGXjwGdYz1XRCpPLGs4/Pa995jTuTOHzJ/Pu1OmcOgdd0BGRoIilGSIZ4KYDbQ3szZmVgvoD+SEH2Bm4XP99gIWhp6/DpxiZg3MrAFwSmibiFSysnop7bFjB1ePHMn0E09kY/36zPjvfzn+rLMSGKEkS9x6Mbl7gZldQfDBngE84e7zzew2INfdc4BhZtYLKADWAReFzl1nZrfzv6VNb3P3dfGKVaS6KquX0kFffcUTF19M1w8/ZFnv3rQbN4529eolMEJJJnOvOtMsZWVleW5ubtkHilRzZU2XsceOHVwzciS33XIL22rXZv7o0Rx7/vlgkZoHJZ2Z2Rx3z4q0T3MxiVQzZY2K7rBgAU8OHsyvP/6YZb17c8CYMRyrld+qpbiOpBaR1BKtvSGjoIAb//535h55JO2WLOH9Z57hgKlTtSxoNaYShEg1Ea294bDPP+fJwYPp/MknPH/OOfhDD3F206YJjlBSjUoQIlVctPENNbdt45ZbbyU3K4sDli3jsueeY/Nzzyk5CKAShEiVFq29odPcuTw5eDCdPvuMZwcMYI8HHmCMJtmTMCpBiFRB0UoNtbZu5fbhw5l99NE0W7mSs6ZOZfszz2gGVilBCUKkiihMCgacT+QurFmzZzOnc2eG33knz5x3HlkLFnB2nz6aZE8iUoIQqQKKT5VRfHRTk9WrefCKK5h1zDE0+PFHer70EreMH8/dDRsqOUiplCBE0lhZE+zV+eUXht9+O0vateOPY8fy6KWXcvT8+ZzXsyd5oOQgUamRWiRNRWuAzigoYPCTT3LriBE0X7GC5888kxv/8Q+2/upX3IsSg8RGCUIkDRUOeCsxpsGdM6ZN464bbqDjwoV8+JvfcO5zz/Fp165ko8Qgu0dVTCJpoLAqaQ+gMXAxJZNDl//+l5ndupHTuzcZO3Zw9pQpdP3gA75XcpByUoIQSXHhDdAOrAW2he1vt3gx/+7bl/8ecwy/+uorLnvkEXrOm8dZZ52Fm6mtQcpNCUIkRZXVAN14zRpGDxvGwg4d6Pnyy/xtxAgOX7yY3152GYtr1lRSkD6V/qUAAA8tSURBVApTG4RICgmfhtuIvM7uXps2cfWoUfz17rups2kT//rDH7h1xAjW7L8/41BpQSqPEoRIiijeK6l4cqj3008Mzc7mqvvvp/mKFUzt04cb//EPvjr4YOqAkoNUOiUIkSQra/GeVnl5XHX//Qx5/HH23biRt088kfOefZZ5v/0t64BWwJ0oOUjlU4IQSaJoYxmO/vhjrr3vPs75z3/YucceTO7fn5HXXMOPRx6phCAJoQQhkgSllRr22LGDM6ZN49r77uO4999nfd263HfttTwwbBg/Zmaqu6oklBKESAIUJoTvgIbABnbtqrrXpk1cOG4cV48axUGLFpHXqhVXjRrFE0OGsGHffWkFSg6ScHFNEGZ2KjAayAD+5e53Fdt/DfAHoABYA1zs7t+G9u0Avggd+p2794pnrCKVrbQeSWvDjtlv1SqueOghLhszhsZr1/Lx0UfTb/Jkppx9Npk1ajAGJQVJnrglCDPLAB4GTgbygdlmluPuC8IOmwtkufsmM7sMuAfoF9q32d07xSs+kXgqq0fSIfPmcdX993P+009Tc/t2cnr14r5rr+X93/6WOmbqkSQpIZ4liC7AYndfCmBmk4HeQFGCcPd3wo6fRTAmSCRtReuR1OiHHxgwaRIXjhtH1pw5bK5dmycuvphRV1/NooMOAtQjSVJLPBNEC2BZ2Ot84NdRjh8CvBr2uraZ5RJUP93l7i9EOsnMhhL8s0bLli0rFLBIeUQb3FZz2zZ6vvwyF4wfT8+XX6bW9u3M7dSJq0aNYsKgQawNreJWB7UxSOqJZ4KwCNsiDQzFzAYBWUC3sM0t3X25mbUFppvZF+6+pMQF3bMJ/rbIysqKeH2RyhStwdkB3MnKzeWC8eMZMGkSjdeuZWXTpjwwbBjjL7iALw4/nJpAXYI/kpao1CCpKZ4JIh84IOx1JrC8+EFmdhLB31s3d99auN3dl4e+LjWzGcCRQIkEIZJIxdsWwhucm3//PYMmTODCcePouHAhW/bckxf69GHchRfy5skns7NGDRxVI0n6iGeCmA20N7M2wPdAf+C88APM7EjgUeBUd18dtr0BsMndt5pZY6ArQQO2SMKFlxj2YNdptvfatIkzp07lwnHjOOmtt9jDnfe7dmXoo4/ybN++rK9fH1BSkPQUtwTh7gVmdgXwOkE31yfcfb6Z3QbkunsOcC+wD/CcmcH/urN2AB41s50Ef5N3Fev9JBJXpbUr7ADqrl/P719/nV45OfTKyaHuhg3ktWrFHcOHM/6CC1hy4IFF11HbgqQzc6861fZZWVmem5ub7DAkDZU1kK1VXh5nTJtGr5wcTpgxg5oFBfzQqBEv9u7N+Asu4L3jjsP3CGbPL0woKjVIOjCzOe6eFWmfRlJLtRQtIawFbOdOusyeTa+cHM6YNo3DvwjGbC7o0IGR11xDTq9ezDrmGHZmZFAzdI11qMFZqhYlCKkWykoIELQnnPTWW/TKyeH0l16i2apVFGRk8N5xx3HNffcx7YwzWNy+PRDUmaqUIFWdEoRUaROBK9m1t1HRc3da5+Vx8ptvcsa0aZz01lvstWUL6+vW5dUePcjp1YvXTj2VHxs23OWaaleQ6kIJQqqUqG0J7rRbsoRuM2dywowZdJs5k5bLgrGc37RuTfbQoeT06sV7xx3H9lq1drmu2hWkOlKCkLQWterInfaLFu2SEDK//x6AVfvtx8xu3bjrhhuYccIJLOzQAex/YzsLB7KpXUGqMyUISXnhSaAlcBrwChFmSXXnV199tUtCaL5iBQArmjVjZrduzDjhBGZ268aXBx+shCBSBiUISTnRSgXfAmNCz+uuX89Rn3xC5zlz6PLxxxz/7rs0W7UKgOX771+UDGaccAJfH3SQEoLIblKCkIQrngDgfx/UpwHjKDmVxd4bN3Lk3Llk5eYWPX719ddF1/y2ZUveOumkooSw+MADd0kI4RoRLFKihCASnRKExFWk6qFICQCC0sFYoPamTfzm00/pPGdOUTLosHAhe4QGdS7LzCQ3K4unzz+f3Kws5nTuzA9NmkR8f5UURMpPCUIqrLQSQaTqobH8r83Adu6kdV4eHRcsoOOCBRwyfz5Hzp1LxwULqLEjmPFoZdOmzD76aJ7t27coGaxq1qzUWJQQRCqPEoSUSyzLaRY+32PHDlrn5XHI/PlFiaDjggV0WLiQOps3Fx3/ffPmfH744bzYuze5WVnkZmWxvHnzUquKCqkLqkh8KEFICaX1GiqtmqgwOdTevJk233zDQV9/vUup4OAvv2SvLVuKrr8sM5MFHTsy9o9/ZEHHjkWPwplPSyy6w66lguLxKCmIxIcSRDW1O9VChb2GcOeXH35g7pIl9F66lHZLltA29LXdkiW0WL7rch/ftmzJgo4dmX7iicw/5BAWdOzIwg4d2FCvXuSVowhGKV+IEoBIKlCCqIJ2twQQXi20ZeNG2uTn0/K774o+/MO/1t2wYZf3+r55c5a0a8ebJ5/MknbtWNq2LYsPPJAFHTuycd99S8RWB/hjWDzFezEpGYikDiWINFD8A//O0PbdLQHss2EDdZYtY0l+Pv3y8zlg2TIy8/OLHgcsW0b99et3ee+ttWqxtG1blrZty7vHH8/Stm1Z0q4dS9q1I691azbXqRM1djUai6QvJYgE2J0P+OL/8Uf6wB9MUE+/DcAdX7eOpqtWcdiqVTQt9ghPAPV+/rlEbCubNiU/M5Ml7doxs1s3lh1wAPmZmeRnZrK0bVu+b9GiaJ2DcMXbCQpfNyp2L0oIIulLCwbFWfE1jCH4r7roA744d/b+5RcarV1b9GiyZk2JD/5mK1fSdNUq9lu9mpoFBSUus71GDVbvtx/ft2hBfmbmLh/8hY/lzZuzbc89d/ue1E4gUnVowaAkugmou2IFRy1evMuHfuGj4bp1JV7vuS1i6mBbzZqsatqUlc2asbx5c+YeeSSrmjaN+PixQYOI//kXKl4CCKdeQyICShBx9x1wc3Y2t/7tb7ts31qrFmsbNSp6fH3QQbu8Dn/80Lgxq5o25af69cscExCL4iUANRSLSCRKEHHWEnjmvPP4oGtX1jVsWPSh/8vee5f7wz5qFRUqAYhI5YhrgjCzUwnmRcsA/uXudxXbvycwHuhM0Nuyn7vnhfbdCAwBdgDD3P31eMYaL3cCQ9u3L1qqEoIP8FqU/gFfXKSeQFB6I7cSgIhUhrglCDPLAB4GTgbygdlmluPuC8IOGwL86O4Hmll/4G6gn5l1BPoDhwDNgbfM7CB33xGveOOl8IO6vL2Yon3gKwmISDzFswTRBVjs7ksBzGwy0BsITxC9gb+Fnv8HeMjMLLR9srtvBb4xs8Wh630Ux3jjZiD6gBeR9FN6N5eKawEsC3udH9oW8Rh3LwDWE3Slj+VcAMxsqJnlmlnumjVrKil0ERGJZ4KI1AJbvGdlacfEcm6w0T3b3bPcPatJKWsCiIjI7otngsgHDgh7nQksL+0YM6sB1COoio/lXBERiaN4JojZQHsza2NmtQganXOKHZND0CUf4BxgugdDu3OA/ma2p5m1AdoDH8cxVhERKSZujdTuXmBmVwCvE3RzfcLd55vZbUCuu+cAjwNPhxqh1xEkEULHPUvQoF0AXJ6OPZhERNKZ5mISEanGos3FFM8qJhERSWNVqgRhZmsIZsRORY2BH5IdRBzp/tKb7i+9VeT+Wrl7xC6gVSpBpDIzyy2tGFcV6P7Sm+4vvcXr/lTFJCIiESlBiIhIREoQiZOd7ADiTPeX3nR/6S0u96c2CBERiUglCBERiUgJQkREIlKCqARmdqqZfWVmi83shgj7jzezT8yswMzOKbbvQjNbFHpcWPzcVFDe+zOzTmb2kZnNN7PPzaxfYiOPTUV+fqH9dc3sezN7KDER754K/n62NLM3zGyhmS0ws9aJijtWFby/e0K/nwvN7IHQejQpJYb7uyb0s/nczN42s1Zh+yr2+eLuelTgQTDP1BKgLcFKop8BHYsd0xo4nGB51XPCtjcEloa+Ngg9b5Dse6rE+zsIaB963hxYAdRP9j1V1v2F7R8NPAM8lOz7qez7A2YAJ4ee7wPUSfY9VeLv57HAB6FrZBAsSHZCsu+pHPfXvfDnAlwG/Dv0vMKfLypBVFzRynnuvg0oXDmviLvnufvnwM5i5/4eeNPd17n7j8CbwKmJCHo3lPv+3P1rd18Uer4cWA2k2qIdFfn5YWadgabAG4kIthzKfX+hpX9ruPuboeM2uvumBMUdq4r8/ByoTfDBuyfB8u+r4h/ybonl/t4J+7nMIlgeASrh80UJouJiXv2uks9NlEqJ0cy6EPwhLqmkuCpLue/PzPYA7gOuj0NclaUiP7+DgJ/M7Hkzm2tm94bWmk8l5b4/d/8IeIegZLsCeN3dF1Z6hBWzu/c3BHi1nOeWoARRcTGvflfJ5yZKhWM0s/2Bp4HB7l7iv/Akq8j9/Ql4xd2XlXlk8lTk/moAxwHXAUcTVHNcVDlhVZpy35+ZHQh0IPiPuwVwopkdX4mxVYaY78/MBgFZwL27e25plCAqriKr36XDynkVitHM6gIvA8PdfVYlx1YZKnJ/vwGuMLM84J/ABWZ2V+WGV2EV/f2cG6reKABeAI6q5PgqqiL3dyYwK1R1tpHgP+9jKjm+iorp/szsJOAmoJe7b92dc6NRgqi4WFbOK83rwClm1sDMGgCnhLalknLfX+j4qcB4d38ujjFWRLnvz90HuntLd29N8F/2eHcv0cskySry+zkbaGBmhe1GJxIs4pVKKnJ/3wHdzKyGmdUEugGpVsVU5v2Z2ZHAowTJYXXYrop/viS7lb4qPIDTgK8J6tdvCm27LfQDg6B4ng/8AqwF5oedezGwOPQYnOx7qcz7AwYB24FPwx6dkn0/lfnzC7vGRaRgL6ZK+P08Gfgc+AJ4CqiV7PupxN/PDIIP1oUEiW9ksu+lnPf3FkHjeuHfWE7YuRX6fNFUGyIiEpGqmEREJCIlCBERiUgJQkREIlKCEBGRiJQgREQkohrJDkCkqjKzRsDboZfNgB3AmtDrLh7MrSOSstTNVSQBzOxvwEZ3/2eyYxGJlaqYREQkIiUIERGJSAlCREQiUoIQEZGIlCBERCQiJQgREYlI3VxFRCQilSBERCQiJQgREYlICUJERCJSghARkYiUIEREJCIlCBERiUgJQkREIvp/uhNX7ouSOgkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#hadron density\n",
    "ax.cla()\n",
    "data = pd.read_csv('./data/hadron_density_T.csv')\n",
    "T_lqcd = data['T'] / 1000\n",
    "rho_lqcd = data['hadron_density']\n",
    "ax.scatter(T_lqcd, rho_lqcd, color='cyan', label='LQCD')\n",
    "ax.plot(T, Nh/0.19732**3, color='red', label='hrg')\n",
    "AddLabelAndLegend(ax, 'T', 'Hadron Density')\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
